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Abstract 

The  standard  recurrence  scheme  does  not  always  yield  accurate  divided  differences  in 
finite  precision  arithmetic.  When  the  function  of  interest  is  known  analytically  and/or  its  values 
are  easily  calculated,  methods  other  than  the  recurrence  scheme  can  be  used.  In  particular,  a 
table  of  divided  differences  can  be  regarded  as  a  function  of  a  special  bidiagonai  matrix.  For¬ 
mulas  and  computational  techniques  suitable  for  computing  matrix  functions  may,  thus,  be 
exploited  for  divided  differences. 

Divided  difference  tables  of  the  exponential  function  are  profitably  treated  as  the 
exponential  of  a  special  matrix.  This  approach  is  good  precisely  when  the  standard  recurrence 
is  bad,  namely  when  the  abscissae  of  the  divided  differences  are  close.  When  the  abscissae  are 
scaled  down  by  powers  of  2,  the  resulting  scaled  divided  difference  table  may  be  squared  to  give 
the  wanted  table.  For  real  abscissae  this  scaling  and  squaring  technique,  in  combination  with 
the  standard  recurrence  where  suitable,  yields  a  hybrid  algorithm  which  permits  computation  of 
any  exponential  divided  difference  to  an  accuracy  dependent  only  on  the  order  of  the 
difference.  For  appropriate  arrangements  of  complex  abscissae,  such  as  conjugate  pairs,  a  simi¬ 
lar  result  is  established.  A  good  way  to  compute  the  exponential  of  a  real  square  matrix  A  is  to 
use  the  Newton  divided  difference  interpolating  polynomial.  Our  algorithm  finds  an  important 
application  in  computing  accurately  the  coefficients  of  this  polynomial.^ 
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ACCURATE  COMPUTATION  OF  DIVIDED  DIFFERENCES 


A  function  of  a  matrix,  f{A ),  may  be  defined  in  terms  of  a  polynomial  which  interpolates 
/  at  A’s  eigenvalues.  One  such  interpolating  polynomial  is  derivable  from  Newton's  divided 
difference  interpolation  formula.  Coefficients  in  this  interpolating  polynomial  are  divided 
differences  of  /  at  the  eigenvalues  of  A.  Thus  /(A )  may  be  represented  in  terms  of  divided 
differences  of  / 

The  opposite  is  also  true,  though  this  is  not  widely  known.  That  is,  divided  differences  of 
/  may  be  represented  in  terms  of  the  function  of  a  special  matrix.  Matrix  functions  and  divided 
differences,  then,  are  profitably  studied  together.  In  particular,  techniques  used  to  compute 
matrix  functions  may  be  exploited  to  study  and  calculate  divided  differences.  The  exploitation 
of  matrix  function  theory  for  the  study  of  divided  differences  is  the  prime  purpose  here.  In  a 
number  of  cases  it  will  lead  to  new  methods  for  accurate  computation  of  divided  differences. 

The  first  chapter  is  a  brief  introduction  to  matrix  functions.  The  interpolating  polynomial 
definition  leads  immediately  to  several  matrix  theoretic  properties  of  /(A),  for  example  A  and 
f(A )  commute.  The  Newton  divided  difference  polynomial  explicitly  shows  the  use  of  divided 
differences  in  defining  f(A).  An  extension  to  a  divided  difference  series  representation  of 
f(A )  is  given  for  hoiomorphic  f 

The  second  chapter  is  a  general  study  of  divided  differences.  §2.1  introduces  a  new  com¬ 
pact  divided  difference  notation  and  lists,  in  this  new  notation,  a  number  of  facts  about  divided 
differences.  For  completeness,  the  following  sections  outline  the  classical  approach  to  the  study 
of  divided  differences  and  the  advantages  of  an  entirely  different  view  of  them  as  functions  of 
their  data  points.  §2.6  establishes  the  matrix  function  formula  for  divided  difference  tables. 
The  remaining  sections  exploit  this  formula  to  develop  series  expansions  for  divided 
differences. 

Chapter  3  is  a  study,  in  detail,  of  divided  differences  of  the  exponential  function  and 
methods  for  computing  them.  The  special  nature  of  /-exp  gives  its  divided  differences  pro¬ 
perties  not  shared  by  those  of  other  functions.  These  properties  are  presented  in  §3.1.  §3.2 
develops  bounds  on  exponential  divided  differences  with  real  data.  These  bounds  show  how 


errors  grow  in  computing  divided  differences  by  the  standard  method.  The  following  sections 
present,  with  error  analyses,  a  Taylor  series  algorithm  and  a  scaling  and  squaring  algorithm  for 
computing  exponential  divided  differences.  The  latter  is  a  direct  consequence  of  representing 
the  divided  difference  table  as  an  exponential  of  a  matrix.  §3.5  then  outlines  a  hybrid  of  those 
two  algorithms  and  shows  how  real  exponential  divided  differences  can  be  computed  with  a 
bounded  relative  error.  Of  prime  importance  is  the  fact  that  the  error  bounds  depend  only  on 
the  order  of  the  difference,  not  the  data.  Finally,  the  remaining  sections  study  complex 
exponential  divided  differences,  with  particular  attention  paid  to  methods  for  computing  divided 
differences  with  data  consisting  of  conjugate  pairs. 
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t.  Matrix  Functions 


1.1  Definitions  and  representations  of  matrix  functions. 

The  extension  to  matrices  of  the  concept  of  function  has  led  to  several  definitions  of 
matrix  function.  Nonetheless,  Rinehart  [19551  has  shown  that  aii  common  definitions  of 
matrix  function  are  equivalent  for  functions  holomorphic  on  a  region  containing  the  eigen¬ 
values  of  the  matrix.  Since  here  we  concentrate  on  holomorphic  functions,  we  are  free  to 
choose  a  definition  which  makes  presentation  easiest.  A  definition  of  a  function  of  a  ma  -tx  in 
terms  of  interpolating  polynomials  has  a  natural  relationship  with  divided  differences.  We 
choose  this  as  our  primary  definition. 

Let  A  be  an  (/;+!)  x  (h-M)  constant  matrix  whose  elements  may  be  complex  numbers. 
We  display  the  eigenvalues  of  the  matrix  A  in  the  sequence 

A, |  =  Ufl . X0,'V| . X. . X, . X /}  m  which  /+l  of  the  eigenvalues  are  distinct 

/ 

and  each  distinct  eigenvalue  occurs  »,+l  times.  /=G,  1 . /.  \4  has  2l<m,  +  1  >  =  w  +  ! 

,«n 

entries.  The  elements  of  A  4  are  just  the  roots  of  A' s  characteristic  polynomial 

x^:xi  -  tx  —  x0)"<)~,(a  —  a '*!  •  ■  •  (x  —  ' .  d.i.i' 

The  definition  we  give  for  /(A)  requires  simply  that  /<X>  be  defined  for  each  X  €  A  , 
when  the  eigenvalues  are  ail  distinct.  To  allow  for  multiple  eigenvalues,  however,  we  lecuire 
that  t  be  defined  on  A  ,  as  follows. 

Definition:  The  function  ./  is  said  to  be  "defined  on  the  characteristic  values  of  .4"  when  /  <  a  ). 

/ ’ ( x  ) . are  defined  for  each  :  «0. ! . i.  For  brevity,  we  denote  this  sequence  o; 

values  by  / t'A  ,). 

For  any  /  satisfying  this  definition,  ./  <  .4 )  is  defined  in  terms  of  an  interpolating  polyno¬ 
mial  for  ./. 


A 


The  polynomial  p  is  an  osculating  interpolation  polynomial  for  /  on  A .  That  is, 
p(A,)-/( A,),  p'(A,)  —/'(A,) . p<"',(A1)— /'"'’(A,)  for  each  /  —  0, 1 . /.  When  the  eigen¬ 

values  are  distinct  this  definition  of  /(A )  becomes  particularly  simple,  as  then  p  is  just  an  ordi¬ 
nary  interpolating  polynomial  for /at  the  elements  of  A^. 

The  rationale  behind  definition  (1.1.2)  is  that  for  two  functions  / and  g,  J'(A)  is  indistin¬ 
guishable  from  gC4)  when  f(\A)  —  g(\A).  The  sequence  of  zeros  of  /(A)-jr(A)  includes 
\A,  the  roots  of  and  xa(A  )  -0  by  the  Cayley- Hamilton  theorem.  The  interpolating 

polynomial  p  has  degree  at  least  />,  since  it  must  satisfy  the  n+l  conditions  given  in  the 
definition/  An  interpolating  polynomial  p  may  be  chosen  to  satisfy  additional  conditions,  but 
the  degree  of  the  polynomial  is  increased.  We  write  p„  for  the  unique  polynomial  of  least 
degree  interpolating /on  A^. 

p„  need  not  be  the  polynomial  of  least  degree  defining  f(A ).  The  characteristic  polyno¬ 
mial  xA  is  an  annihilating  polynomial  for  A  because  xa  (A )  ”0.  However  for  some  matrices  A . 
there  are  polynomials  of  smaller  degree  which  are  also  annihilating  polynomials.  The  minimal 
polynomial  p.A  is  the  non-trivial  annihilating  polynomial  for  A  of  least  degree.  If  /14(A)  has 
degree  m+l,  m  <  n,  it  is  possible  to  define  f(A )  in  terms  of  a  m  degree  polynomial  pm  which 
interpolates  /  at  the  m+1  roots  of  fiA.  Gantmacher  [19591  uses  this  slightly  more  general 
approach  in  his  definition  of  f(A ).  The  roots  of  nA(k)  are  eigenvalues  of  A.  For  m  <  n  fewer 
derivatives  of  /  need  be  specified,  however  nA  and  the  multiplicities  of  its  roots  may  be  difficult 
and  costly  to  obtain.  Thus  we  shall  not  try  to  form  f(A)—pm(A)  for  the  smallest  possible 
degree  m.  p„  can  have  significantly  higher  degree  than  p„,  see  Fig.  1.2.1,  but  here  we  achieve 
greater  simplicity  in  that  less  need  be  known  about  the  matrix  A. 

|A  polynomial  of  degree  k  can  interpolate  at,  at  most,  k+ 1  points.  In  general  k+ 1 
points  uniquely  determine  a  polynomial  of  degree  ic,  higher  degree  polynomials  are  not 
uniquely  determined. 


0  12]  M^(X)  -  (X  —  1 ) (X  —  2) 

A  -  0  1  0 

-1  1  3J  -  (X-  1)2(X  — 2) 

Fig.  1.1.1:  Degree  of  nA  may  be  less  than  degree  of  *,«. 


The  polynomial  representation  of  f(A)  leads  to  several  elementary,  but  very  useful, 
consequences. 


Similarity  transformations.  For  any  (n+1)  x  (n+1)  nonsingular  matrix  P, 

APAP-')  -  P  f(A  )•/»-'. 


(1.1.3) 


In  theory  this  permits  performing  all  computations  to  form  f{A )  on  the  simplest  matrix  similar 
to  A,  e.g.  /4’s  Jordan  canonical  form.  In  practice,  however,  the  transformation  matrix  P  may  be 
difficult  to  compute  accurately1,  or  may  be  nearly  singular.  Some  less  simple  form  may  be 
required.  The  triangular  Schur  form  T,  which  is  unitarily  similar  to  A,  eliminates  the  above 
objections.*  However,  f(T)  is  not  always  simple  to  compute  with  accuracy. 


Commutativity. 

A  f(A)  -AA)-A 

(1.1.4) 

Parlett  [1976]  has  presented  a  very  fast  method  for  computing  functions  of  upper  triangular 
matrices  T  based  on  this  property.  In  brief,  the  diagonal  of  f(T),  which  is  also  upper  triangu¬ 
lar,  is  computed  directly; 

/(r),,-/(7V) 

for  each  i-0,1 . n.  Then  successively  by  diagonals  towards  the  upper  right,  the  general 

recurrence  is 

fKagstrdm  and  Ruhe  [1976]  present  an  algorithm  for  computing  the  Jordan  form, 
while  Golub  and  Wilkinson  [1976]  discuss  limitations  on  computing  it  accurately. 

♦Wilkinson  [1965]  presents  a  detailed  analysis  of  the  QR  algorithm  which  reduces  A  to  . 

T  by  a  sequence  of  unitary  similarity  transformations;  the  algorithm  is  implemented  in 
the  EISPACK  [Smith,  1974]  collection  of  computer  subroutines. 
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AT )tJ  -  |rfV(n,+rU  -  r„_*-/(n/.*J))/(ru-  tj  u.i.5) 

t-0 

where  /  <  j  <  n.  T  may  be  A's  Schur  form;  this  recurrence  may  be  used  to  form  f(A )  by 
(1.1.3). 

When  /  is  symmetric  in  the  real  axis,  that  is  /({)  -/({),  polynomials  interpolating  /  have 
real  coefficients.  We  denote  the  conjugate  transpose  of  A,  A T,  by  A”. 

Conjugate  transpose.  When  /  is  symmetric  in  the  real  axis, 

/CO  -/CO*.  (1.1.6) 


i 


\ 

> 


t 


Expression  (1.1.6)  shows  that  conjugate  symmetries  in  A  are  inherited  by  /(A). 


Formula  (1.1.3)  shows  that  f(A ),  defined  as  in  (1.1.2),  may  always  be  computed  from 
/4’s  Jordan  canonical  form.  Conversely,  we  may  wish  to  define  f(A )  from  the  Jordan  form  by 
way  of  (1.1.3).  This  latter  definition  is  more  general  than  our  polynomial  definition,  as  the  fol¬ 
lowing  shows. 

The  2x2  identity  matrix  has,  among  others,  the  square  roots 


1  0 
0  I 


and 


1  0 
0  -1 


The  former  root  is  representable  by  either  definition;  the  latter  is  obtained  by  separately 
defining  VT-1  and  VT--1  on  each  Jordan  block.  The  function  is  permitted  to  be  mul¬ 
tivalued,  but  only  on  separate  Jordan  blocks.  The  polynomial  definition  does  not  allow  this, 
since  polynomials  are  never  multivalued. 

Even  the  Jordan  form  definition  of  f(A )  is  not  the  most  general  possible.  For  example  a 
square  root  of 


0  0 
0  0 


is 


0  1 
0  0  • 


E.  Cartan  proposed  a  contour  integral  definition  which  applies  to  holomorphic  functions  / 
[Rinehart,  1955]. 


§1.1 

Cartan  definition.  If  /  is  holomorphic  inside  and  on  a  simple  closed  contour  C  enclosing 
A^,  then 

/(A)  s  (1.1.7) 


Additional  representations  of  f(A )  are  derivable  from  those  just  mentioned.  Gantmacher 
[1959]  and  Rinehart  [1955]  discuss  f(A)  in  further  detail.  In  the  next  section  we  present  a 
particular  polynomial  representation  for  f(A )  and  discuss  related  series  representations. 


1.2  The  Newton  polynomial  of  f(A),  and  series  representations. 

Because  we  are  free  to  choose  any  polynomial  interpolating  /  on  A A,  definition  (1.1.2) 
allows  many  representations  of  f(A).  There  is,  however,  a  unique  interpolating  polynomial  p„ 
of  least  degree,  though  even  this  may  be  arranged  in  many  ways/  One  arrangement  of  p„, 
which  clearly  illustrates  the  use  of  divided  differences  for  defining  matrix  functions,  is  based 
upon  Newton’s  divided  difference  formula  for  the  interpolating  polynomial,  namely 

P„(k)  -  lAo‘/  n(X-X() .  (1.2.1) 

4*0  7— 0 

The  coefficient  A<j/  is  the  k- th  order  divided  difference  of  /  defined  on  the  abscissae 
A0.  A, . A*.  This  compact  divided  difference  notation  is  further  explained  in  §2.1. 

The  first  few  terms  of  the  interpolating  polynomial  (1.2.1),  which  we  call  a  Newton  poly¬ 
nomial,  are 

/(Xo)  +  Ao7(A-A0)  +  A^/-(X-Xo)(X-X,)  +  A(J/(X-X0)(X-X,)(X-X2)  +  ■  •  • 

Because  p„( \A)-f(\A)  where  A4-{X0.X| . A„),  the  eigenvalues  in  A*  having  been 

renumbered,  f(A)  has  the  following  representation. 

Newton  polynomial  of  f(A).  When  /is  defined  on  the  characteristic  values  of  A, 

f(A)  -  3>«?/-n(-4-A,/).  d.2.2) 

4—0  7—0 

A  A  is  the  sequence  of  abscissae  for  the  divided  differences. 

In  §2.1  we  shall  see  that  the  conditions  on  /  necessary  to  define  all  the  divided  difference 

coefficients  A <f/,  4—0, 1 . n,  are  exactly  those  required  to  assure  the  existence  of  some 

interpolating  polynomial  p„.  Thus  when  p„  exists,  it  may  be  arranged  as  a  Newton  polynomial; 
so  (1.2.2)  is  equivalent  to  definition  (1.1.2). 

ft 

fFor  example  Lagrange's  interpolating  formula  p„(k)  -  X/*(A)-/(A*),  where  for  each 

*- 0 

k  lk(k)  =  r|(X-X^)/J3(XA -Xy),  is  one  of  the  simplest.  Here  /*(X,) -0  when  k  /, 
y-o  j  *k 

and  /*(A*)  - 1. 
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1 

1  1 

1 

1  1 

1 

B  - 

1  1 

1 

1 

2 

2 

M,*<X)  -  (X-lMx-2)  MaOO  -  (X  —  1)4(X  —  2) 

X^(X)  —  (X  —  1)4(X  — 2)  X*00-M*(X) 

PM )-/(!)/  +  [f(2)-Al)](A  -/)  pJB )  -  S  (g-/)*  +  A0V(g-/)4 

1-0  *  • 

3 

Zi-~L(^-/)‘  +  Ao4/(><-/)4  p„(B)-pJB) 

*- o 

Fig.  1.2.1:  p„,  depends  on  the  eigenspaces  of  the  matrix. 

The  Newton  polynomial  representation  of  /(.4 )  requires  no  more  of  /  than  that  it  have 
enough  derivatives  to  define  /(A A).  Our  interest  here,  however,  concerns  functions  /holo- 
morphic  on  a  region  containing  \A.  In  such  cases  there  is  a  natural  extension  of  the  Newton 
polynomial  to  a  series.  Such  a  series  may  be  viewed  as  an  interpolating  polynomial  of  infinite 
order. 

This  extension  derives  from  a  Newton  divided  difference  series, 

/(X)  -  £A<f/-n(X-My) .  (1.2.3) 

*-0  7-0 

where  the  divided  differences  of  /  are  defined  on  a  sequence  of  expansion  points 
M  =  !  which  lies  in  the  domain  of  holomorphy  of  f.  Because  (1.2.3)  may  be 

unfamiliar.  Appendix  A.l  presents  an  elementary  proof  demonstrating  its  convergence.  Appen¬ 
dix  A.2  establishes  the  following  representation  of  f(A).  Gantmacher  (1959]  establishes  more 
general  series  representations  for  f(A),  and  Gel’fond  [1971]  discusses  more  complicated 
divided  difference  expansions.  These  more  extensive  results  are  not  needed  here. 
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Newton  series  representation  of  f(A).  Let  /  have  the  Newton  divided  difference  expan¬ 
sion  (1.2.3)  on  an  open  disk  containing  A,*.  Then 

f(A)  -  £A0‘y  nC4 (1.2.4) 
*« 0  /«0 


When  the  first  /j+1  elements  of  M  comprise  i.e.  . the  Newton  expansion 

of  /(A)  (1.2.4)  terminates  after  the  n-th  term  and  is  just  the  Newton  polynomial  (1.2.2).  This 

H 

is  the  Cayley-Hamilton  theorem,  nu-M)-^(^)-o. 

7-0 

When  M- (m.M.M.— )  consists  of  one  point,  then  each  (§2.1).  The 

Newton  expansion  (1.2.3)  is,  then,  just  a  Taylor  series;  the  representation  (1.2.4)  reduces  to  a 
Taylor  series  for  f(A). 

Taylor  series  representation  of  f(A).  Let  /  have  a  Taylor  series  on  an  open  disk  about  n 
containing  A^.  Then 

f(A)  -  £  (A  -*/)*.  (1.2.5) 

7-0  *■ 


The  above  shows  that  f(A)  is  representable  in  terms  of  f's  divided  differences.  In  the 
next  chapter  we  reverse  this  situation.  Divided  differences  of  /  are  expressed  in  terms  of  a 
function  of  a  special  matrix.  Hence  everything  said  here  concerning  f(A )  applies  to  divided 
differences  of  f,  and  techniques  suitable  for  computing  f(A )  may  be  applied  to  compute  them. 
In  turn,  these  differences  may  be  used  to  compute  f(A )  by  the  Newton  polynomial. 


n 


Si  1 


2.  Divided  Differences 


2.1  Definitions  and  properties  of  divided  differences. 

Divided  differences  were  studied  extensively  in  classical  precomputer  numerical  analysis 
as  part  of  a  finite  difference  calculus.  They  primarily  saw  use  in  tabulation  of  tables  of  function 
values.  A  quite  different  purpose  is  envisioned  here;  however,  much  of  the  classical  theory  is 
still  relevant.  Before  proceeding  to  develop  formulas  for  the  calculation  of  divided  differences, 
we  present  a  few  well-known  definitions  and  their  consequences.  Our  notation  is  somewhat 
different  from  that  of  other  authors,  but  it  is  felt  to  be  an  improvement.  Once  understood,  it 
will  cause  no  confusion  to  those  already  familiar  with  divided  differences. 

Most  common  notations  for  divided  differences  are  cumbersome.  For  clarity  we  begin 
with  such  a  notation,  but  later  reduce  it  to  more  compact  form  by  suppressing  unneeded  infor¬ 
mation.  Let  /  be  a  function  of  a  single  variable  £  and  be  defined,  at  least,  on  a  sequence 

Z-(fo.{| . £„.... )  of  distinct  complex  numbers.  Z  is  called  the  sequence  of  abscissae,  or 

sometimes  the  sequence  of  data  points  or  nodes.  The  0-th  divided  difference  of  /at  £o  is 

(A0/)  (Co)  *  /(Co)  • 


The  first  divided  difference  of  /at  £0  is  a  function  of  the  two  variables  (abscissae)  £0  and  £i;  it 
is  formed  from  the  0-th  divided  difference  by  the  familiar  formula 


(A'/)(£0.<i) 


U°/)(£i)-(A°/)(£o)  /<£i)-/(£q) 

Ci— Co  "  Ci  —  Co 


The  fc-th  order  divided  difference  of  /  at  £0  is,  then,  a  function  of  the  *+l  abscissae 
£o,£i . £*,  and  is  defined  iteratively  from  fc-l-st  order  divided  differences. 


A  first  definition  of  divided  differences.  When  /  is  defined  on  Z,  each  Ar-th  order  divided 
difference  of  /at  £,,  j  -0, 1 . n-k,  is 


UA/)(£,.£y+l . «,+*> 


(A  *->/)(£,+, . -  (A^-'/HC/ . C/+*-i)  (2  n 


(A*/H{>.{,+ 1.  •  •  •  .{>+*)  has  no  dependence  on  abscissae  with  indices  <)  or  >y'+Ac,  and  so 
no  generality  is  lost  when  considering  just  (AYHCo^i . («)- 
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Divided  differences  are  very  special  functions  of  the  data  points  in  Z.  Not  only  does  the 
number  of  data  points  used  increase  with  the  order  of  the  difference,  but  the  divided  difference 
is  symmetric  in  its  arguments.  This  is  obvious  in  the  equivalent  representation  of  the  divided 
difference  in  terms  of  determinants  [Milne-Thomson,  19331. 


/(Co) 

/(c.) 

f(U 

Co 

C." 

r  " 

Vtt 

(AY) (Co . C«)  “ 

cr1 

cr1 

r  »-> 

WM 

cr1 

cr‘ 

. 

i 

i 

l 

i 

i 

1 

The  abscissae  may  be  arranged  in  any  order  without  changing  the  value  of 

(A"/)  (Co.  Ci . C„). 


When  /  is  symmetric  in  the  real  axis,  i.e.  /(C)  “/(C),  (2.1.2)  leads  to  a  conjugate  symmetry. 

For  odd  values  of  n,  (A"/)(Co.Ci . C»)  is  real  whenever  C2/+i“C2m  '“0,1 . (n-\)/2. 

And  for  n  even,  (!"/) (Co- Ci . („)  is  the  conjugate  of  (A"/)(Co.Ci.  •  •  •  -C»)  when  each 

C2,>i-<m  '  “0,1 . (n—2)/2. 

The  defect  in  definition  (2.1.1)  is  that  data  points  must  be  distinct.  However  when  /  is 
differentiable,  (2.1.1)  may  still  be  defined  even  for  confluent  (i.e.  equal)  abscissae.  In  particu¬ 
lar  when  Z-ICo. Co . Coh  (2.1.1)  is  defined  when  /‘"’(Co)  exists.  For  confluent  abscissae 

the  divided  difference  reduces  to 

f<nUrn) 

(A "/)  (Co.  Co . Co) - (2.1.4) 

Since  the  data  points  may  be  arranged  in  any  order  without  changing  the  value  of  the  divided 
difference,  (2.1.1)  is  defined  when  (2.1.4)  is  used  when  confluent  abscissae  occur.  The  require¬ 
ment  that  the  abscissae  be  distinct  may  be  removed. 

Definition:  Let  Zs(Co. . . . . . ii . C/.-l  be  a  sequence  of  abscissae 

gust  a  renumbering  of  the  previous  Z)  where  each  („  '“0,1 . /,  appears  *,+1  times, 

£(n,+l)  -  fl-t-1.  The  function  /  is  "defined  on  the  sequence  of  abscissae  Z"  when  /((,), 
-o 
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./'({,) . /f  ,  ,({,)  are  defined  for  each  /  — 0. 1 . /.  This  sequence  of  values  is  denoted 

/( Z). 

Before  rewriting  definition  (2.1.1)  in  more  generality,  we  introduce  a  compact  notation. 
The  sequence  of  data  points  Z —  (£<>,  £|.  .  .  .  is  given  and,  usually,  in  a  fixed  order. 

Hence  reference  to  Z  may  be  suppressed.  Thus  we  define 


1,7  3  U7>(i/.{/+ . £,>*) . 


(2.1.5) 


The  subscript  j  is  understood  to  mean  that  we  locate  the  abscissa  tabled  £,  and  use  it  and  the 
next  k  abscissae  in  the  sequence,  in  the  event  that  the  particular  sequence  Z  must  be 
emphasized,  A  will  be  written  for  A  ,f.‘r 


Standard  iterative  divided  difference  scheme.  When  /  is  defined  on  the  sequence  of 
abscissae  Z, 


(2.1.6) 


for  each  *-l,2 . nan  d  0. 1 . ff— it,  where  A  °/ =/((,). 


This  definition  of  divided  differences  and  our  earlier  definition  of  matrix  functions  in  §1.1 
are  consistent.  Indeed  when  Z  -  A ,  "defined  on  the  sequence  of  abscissae  Z"  and  "defined  on 
the  characteristic  values  of  A "  are  the  same.  We  shall  see  later  in  §2.6  that  this  similarity  in 
definitions  is  no  coincidence. 


Divided  differences  have  many  useful  representations  and  properties.  We  list  several  of 
these  here. 

Divided  difference  tables.  Divided  differences  are  most  conveniently  displayed  in  tables.  Trad¬ 
itionally,  tables  are  arranged  as  in  Fig.  2.1.1.  Each  divided  difference  is  computed  from  its  two 
immediate  neighbors  in  the  column  to  its  left.  For  our  purposes  it  is  most  helpful  to  arrange 

tMilne-Thomson  (1933]  writes  Ao/ as  [{0,£i . £„],  suppressing  the  function;  Davis 

(1973]  uses  /"'({o.ii . {„);  and  Kahan  and  Farkas  (1963]  use  A/({0,s i . £J, 

which  suggested  the  notation  used  here.  Gabel  [1968]  also  uses  a  similar  notation. 

This  compact  notation  is  used  in  McCurdy  [1978],  from  which  much  of  this  introducto¬ 
ry  section  is  taken. 


i 

•  > 

y 

l  a. 
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Co 

/(Co) 

Ao1/ 

C, 

/((.) 

A  If 

A,1/ 

A  If 

Ci 

nci) 

A,2/ 

AoV 

A,V 

C.t 

/(£j> 

A,1/ 

C4 

AC  4> 

Fig.  2.1.1:  Standard  divided  difference  table. 


the  table  as  an  upper  triangular  matrix,  for  example 


A/  a 


/(Co)  Ao1/ 

A<£/  ■ 

Ao"/ 

/(Ci) 

A,’/  • 

•  AT1/ 

/(c2)  • 

•  Ar2/ 

/<{«> 


(2.1.7) 


The  symbol  A/,  without  the  superscript,  is  used  here  to  represent  a  matrix,  not  a  scalar.  Ele¬ 
ments  of  the  matrix  depend  on  their  immediate  neighbors  in  the  diagonal  to  the  left.  This 
leads  to  a  'pattern  of  dependence”  in  which  A  ff  is  independent  of  all  table  entries  in  rows 
before  the  >th  and  columns  after  the  j+k-\h.  A  */  depends  only  upon  the  block  of  the  table 
matrix  between  it  and  the  main  diagonal.  Such  patterns  of  dependence  are  characteristic  of  tri¬ 
angular  matrices. 

Linearity.  For  constants  a  and  £, 


&A'(af+fig)  -  a-A<j If  +  fl-Lfig  . 


(2.1.8) 


Translation  invariance.  For  Z+a  =  {£o+a,£i+a,  •  •  .{«+<*,... )  and /a({)  =/((+a). 


A{J/o  -  L(W- 


(2.1.9) 


/■(Ci) -/.(Co) 
5i  —  Co 


/(Ci+«)—  /<Co+«)  .  ,  , 

(Ci+a)-({0+«) 


A 


For  example, 
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/({,) 


A,1/ 

-  A,2/  - 

A,3/ 

-  A,4/ 

r 

1 

I 

I 

/(«<♦!) 

-  A, Vi/  - 

A,Vt/ 

-  A,+,/ 

1 

t 

t 

/<C,+i>  - 

A, Vt/ 

-  A, Vi/ 

t 

1 

/({,+>) 

A  ,Vi/ 

? 

/({,+4) 

Fig.  2.1.2:  Pattern  of  dependence  in  a  divided  difference  table. 


Scaling  invariance.  For  some  r^O  let  tZ  =  {tJo,  r£i,  .  .  .  .rf,,,...)  and  /T(£)  =/(t£).  Then 


(2.1.10) 


For  example. 


fAU-fAlo) 

Ct-Co 


/(r{,)-/(TCo) 


1.000  1.718  1.476 

2.718  4.671 

7.389 


.8455  .3632 

4.013  2.298 

12.70  10.91 

20.09  34.51 

54.60 


Fig.  2.1.3:  Divided  difference  table  for/-exp,  with  Z-{0,1,2,3,4}. 


Mean  value  representation.  When  the  abscissae  are  real. 


A0"/ 


/"»({) 

nl 


min  t, 

OiJ^n  1 


<  £  <  max  £,, 
0<< 


(2.1.11) 


for  any  /  having  n  continuous  derivatives  in  the  interval  containing  the  data  points.  This  has 
no  equivalent  for  complex  abscissae.  For  example  when  /—exp  and  £o-f  and  £,  -£+2wi. 


Ao'exp 


et+2  *'—et 

((  +  2iri)-S 


0  ^  e( 


for  any  finite  £. 


rr 


Integral  representation.  Another  representation  for  An'/;  when  /  has  a  bounded  n-th  order 
derivative  on  a  closed  convex  domain  containing  Z,  is  [Gei'fond.  1971] 


i  Ti  T,,-\ 


Atf/-  ff  J/',ICo+«|-Co)T.+  •  •  ■  ■  •  rftjrfr,.  (2.1.12) 


u  o  n 


Contour  integral  representation.  When  /  is  hoiomorphic  inside  and  on  a  simple  closed  con¬ 


tour  C  enclosing  Z,  [Gei’fond,  1971] 


/(at)  dcu _ 

{0)(ti»-£|)  •  •  •  (w-<„) 


(2.1.13) 


Bound.  If  /has  a  bounded  n- th  derivative  on  a  closed  convex  domain  Cl  containing  Z,  then 
[Gei’fond,  1971] 


|Ao/|  <  -Jrmajt  |/"’({)  I  • 
an 


(2.1.14) 


This  is  an  immediate  consequence  of  (2.1.12). 

In  later  sections  we  present  a  new  way  of  looking  at  divided  differences  and  develop  addi¬ 
tional  ways  to  express  them. 


1 
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2.2  A  traditional  attitude  towards  divided  differences:  tables  and  interpolation. 

Before  beginning  our  study  of  divided  differences,  we  present  a  short  discussion  of  the 
traditional  attitude  towards  them  in  order  that  contrasts  may  be  made  with  our  approach. 

Divided  differences  are  often  encountered  as  an  adjunct  of  the  subject  of  ordinary 
differences  in  interpolation  and  table  making.  Their  treatment  in  the  literature  [e.g.  Milne- 
Thomson,  1933,  and  Miller,  1950]  is  patterned  on  that  for  ordinary  differences.  The  arrange¬ 
ment  of  the  divided  difference  table  (Fig.  2.1.1)  is  one  example.  Others  are  divided  difference 
interpolation  formulas  which  resemble  formulas  for  ordinary  differences. f  indeed,  our  borrowed 
notation  A  for  the  divided  difference  operator  is  a  modification  of  A  for  the  foreward  difference 
operator. 


Fig.  2.2.1:  Error  growth  pattern  in  table  of  ordinary  differences. 


The  lack  of  interest  in  divided  differences  shown  by  some  authors  [e.g.  Ralston,  19651  is 

explicable  when  we  recall  their  use  in  interpolation  and  the  available  means  of  computation.  In 

tDivided  difference  interpolation  formulas  are  derived  from  the  Newton  divided 
difference  interpolating  polynomial  (1.2.1)  in  the  same  way  that  the  Stirling  and  Bessel 
formulas  are  derived  from  Newton’s  foreward  difference  scheme.  For  example,  averag¬ 
ing  polynomials  using  the  two  sequences  of  abscissae  {-:••••)  and 

{(o>  C-i- {t>  t-2'  ii'  •••)  yields  the  divided  difference  generalization  of  Stirling’s  formula: 

/*<()  -  m  0>  +  j|(A '  miflO  ♦  «A  '/7<to-  -  W  +  (Aj/)({0.{|.{-|)({-{0><{--^7=1-> 


1  idU  ■ 
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hand  or  desk-top  calculation,  we  desire  few  and  simple  computations.  Ordinary  difference 
interpolation  formulas  may  be  scaled  to  minimize  divisions  and  complicated  fractions  (worse 

than  say  y,  y,  -jy ,  etc.).  This  is  one  reason  mathematical  tables  with  unevenly  spaced  data 

are  seldom  encountered.  Comrie  (1959)  remarks  that  ’computers  try  to  avoid  tabulations  at 
unequal  intervals  and  divided  differences  ...  .’ 


i 

A°exp 

A  'exp 

A2exp 

AJexp 

A4exp 

A5exp 

A6exp 

0.00 

1.000 

0.284 

0.25 

1.284 

0.365 

0.081 

0.122 

0.50 

1.649 

0.568 

0.203 

-0.270 

-0.392 

1.000 

0.75 

2.217 

0.501 

-0.067 

0.338 

0.608 

-0.997 

-1.997 

1.00 

2.718 

0.772 

0.271 

-0.051 

-0.389 

1.25 

3.490 

0.992 

0.220 

1.50 

4.482 

0.00 

0.000 

0.000 

0.25 

0.000 

0.000 

0.000 

0.100 

0.50 

0.000 

0.100 

0.100 

-0.300 

-0.400 

1.000 

0.75 

0.100 

-0.100 

-0.200 

0.300 

0.600 

-1.000 

-2.000 

1.00 

0.000 

0.000 

0.100 

-0.100 

-0.400 

1.25 

0.000 

0.000 

0.000 

1.50 

0.000 

[  ?  Fig.  2.2.2:  Example  of  error  propagation  in  ordinary  differences  (see  Fig.  2.2.5). 

'  I 

The  study  of  error  behavior  in  divided  difference  computations  also  shows  the  domination 
of  ordinary  difference  theory.  For  example  when  an  error  f  occurs  only  in  the  /(0)  entry,  the 
familiar  error  growth  pattern  (Fig.  2.2.1)  reveals  itself  in  a  table  of  ordinary  differences.  The 
coefficient  (exclusive  of  sign)  of  the  (»,  Ay)  entry  is  the  binomial  coefficient 


J 

j/2  -  n 


,  u 
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where  n  is  assigned  half  values  for  entries  where  j  is  odd. 

Miller  [1950]  writes,  "It  is  also  proposed  to  give  error  patterns,  such  as  that  in  [Fig.  2.2.1 1, 
for  tables  of  divided  differences,  for  use  with  tables  having  certain  common  arrangements  at 
unequal  intervals,  for  example,  with  a  table  having  arguments 


1 

3’ 


2' 


1.  I  i.  ,1  ,1 

3  4  4'  3 


N 


Such  an  error  pattern  might  resemble  Fig.  2.2.3. 


Fig.  2.2.3:  Error  growth  pattern  in  table  of  divided  differences. 


Each  (So-C*)l-1-  The  error  growth  pattern  reduces  to  that  of  Fig.  2.2.1  when  the 

k-t 

k*  0 

data  points  are  evenly  spaced  with  unit  separation  and  each  entry  eu  is  multiplied  by  j\. 

Because  w-th  order  differences  of  polynomials  of  degree  «-l  are  zero  (see  §2.7),  one 
expects  high  order  differences  of  a  function  to  be  small  when  it  is  well-approximated  by  a  poly¬ 
nomial.  When  high  order  differences  begin  resembling  the  alternating  sign  binomial  pattern  of 
Fig.  2.2.1,  an  error  in  the  tabulated  function  values  is  suspected  [Mitler,1950l.  Multiple  tabula¬ 
tion  errors  lead  to  more  complicated  patterns,  and  round-off  errors  in  the  difference  computa¬ 
tions  may  further  obscure  any  pattern.  Statistical  methods  have  been  suggested  for  spotting 
aberrations  in  tables  [Blanch,  1954], 
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Fig.  2.2.4:  Example  of  error  propagation  in  divided  differences  (see  Fig.  2.2.5) 


A  situation  where  divided  differences,  rather  than  ordinary  differences,  are  usefully 
employed  in  interpolation  is  presented  by  Salzer  (19471.  Bessel  functions  ),(£).  etc , 

are  commonly  tabulated  for  integral  values  of  v ,  as  well  as  v—  ±1/4,  ±1/3,  ±1/2,  ±2/3,  and 
±3/4.*  For  £  fixed,  divided  differences  are  used  to  interpolate  for  any  v,  -1  $  1,  or  ;o 

check  entries  in  a  table. 

Very  high  (say  greater  than  10-th)  order  differences,  ordinary  or  divided,  are  seldom  of 
practical  interest  in  interpolation  problems.  The  reason  is  that  when  the  function  is  tabulated  to 
a  fixed  number  of  digits,  adjacent  table  entries  often  have  several  initial  digits  in  common. 


fSee  for  example  tables  of  the  National  Bureau  of  Standards  (19481 


J 


_ j*£  *tjMa fee 
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Differencing  such  entries  leads  to  differences  containing  fewer  significant  digits.  After  a  few 
steps  no  correct  digits  may  remain.  The  tables  in  Fig.  2.2. 5  illustrate  this  phenomenon.  Less 
correct  information  (significant  digits)  is  retained  after  every  differencing  step,  in  interpolation, 
one  is  interested  only  in  making  a  small  correction  in  the  last  digits  of  already  tabulated  values. 
The  information  remaining  in  the  first  few  differences  is  adequate  for  this  task.  However  when 
accurate  high  order  differences  are  the  objects  of  interest,  this  loss  of  information,  coupled  with 
magnification  of  any  previously  introduced  errors  when  we  divide  by  a  small  number,  is  a  disas¬ 
ter.  We  must  consider  other  methods  for  computing  divided  differences.  The  approach  neces¬ 
sary  to  develop  such  methods  forsakes  the  idea  of  interpolation  between  table  entries  and 
emphasizes  the  underlying  function. 

example:  We  use  the  Newton  divided  difference  formula  and  the  four  figure  divided 
differences  in  the  second  table  of  Fig.  2.2.5  to  interpolate  for  exp(0.30). 

exp(0.30)  -  1.000+  1.136  x  (0.30-0.00)  +  0.648  x  (0.30-0.00X0.30-0.25) 

+  0.2347  x  (0.30 -0.00)  (0.30 -0.25)  (0.30 -0.50)  +  ••• 

-  1.349816 

This  result  correctly  interpolates  to  four  figures  exp(0.30)  •«  1.350.  The  errors  in  the 
divided  differences  do  not  affect  the  most  significant  digits  that  we  want. 


m'jj' 


5 

0 

00 

.25 

.50 

.75 

lull 


Auexp  A'exp  A2exp 


284  |  0.081 

365 


0.772  0.220 

0.992 


Ordinary  differences  of  the  exponential 


A’exp  A4exp  A5exp 


8  0.000 

8  0.003 

0.011 


Divided  differences  using  4  digits 


A°exp 


1.872  1.064 

2.404  1.368 

3.088  1.760 

3.968 


5.227E-1 


A4exp 

A5exp 

Aftexp 

8.530E-2 

8.530E-2 

1.174E-1 

0.000 

2.568E-2 

1.712E-2 

Correct  value  of  divided  differences  to  4  digits 


A  “exp 


A  *’exp 

AJexp 

A4exp 

6.454E-1 

8.287E-1 

1.064 

1.366 

1.754 

2.444E-1 

3.138E-1 

4.029E-1 

5.174E-1 

6.942E-2 

8.913E-2 

1.144E-1 

Fig.  2.2.5:  Example  of  loss  of  accuracy  in  computing  differences. 
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2.3  An  analytic  approach  to  divided  differences. 

Up  to  this  point,  divided  differences  are  seen  as  a  tool  for  interpolating  in  mathematical 
tables.  We  assume  only  that  we  are  given  a  table  of  numbers,  presumably  representing  the 
values  of  some  function  at  certain  arguments.  No  reference  to  a  particular  function  or  expres¬ 
sion  is  required. 

In  contrast,  we  now  consider  how  divided  differences  depend  on  the  function  / and  make 
full  use  of  the  theoretical  tools  presented  in  §2.1.  We  treat  the  divided  difference  itself  as  a 
function:  hence,  we  always  assume  we  can  evaluate  /  and  its  derivatives  at  any  valid  abscissa.* 
A  discussion  of  tables  and  interpolation  is  no  longer  relevant;  neither  is  a  limitation  to  common 
arguments,  as  suggested  by  Miller  [1950]  and  Salzer  [1947].  Indeed,  complex  as  well  as  real 
data  points  are  possible.  Further,  we  are  interested  in  floating-point  computation  on  a  com¬ 
puter;  the  desire  to  avoid  divisions,  complicated  numbers  (many  digits)  and  fractions  is  less 
important.  Finally  ,  we  consider  divided  differences  of  any  order. 

example:  The  power  of  the  analytic  approach  can  be  illustrated  as  follows.  We  wish  to  evalu¬ 
ate  A ’exp  at  the  abscissae  £O“0  and  IO“20  on  a  pocket  calculator  that  can  hold 
only  a  ten  digit  number.  In  the  calculator  the  number  1  + 10“20  would  be 
1.000000000;  the  10“ 20  is  chopped  off.  Hence  we  compute 


Ao'exp 


exp(10~2°)  -cxp(O) 
10“2°-0 


1-1 

10“ 20 


0 


asexpOO  20) 


A,'exp  - 


•  1  to  ten  digits. 

exp(g|)  -exp(gp) 
ft -Co 


Alternatively  we  may  write 

sinh[(Cl-C0)/2] 


exp[(4i  +  {„)/2]- 


<{i-«o)/2 


and  then 

.i  ,n  c ..  ,„-!()<  sinh<0.5  x  I0~;u) 

A,]exp  -  exp(0.5  x  10  20) - — --ps - 

to  ten  digits,  if  we  can  evaluate  sinh  accurately. 


tin  a  computer,  abscissae  must  be  representable  in  the  machine;  the  value  of  /  is 
rounded  (or  chopped)  at  full  machine  precision. 


Divided  differences  of  g  using  4  digits 


£ 

A0* 

A1* 

A-V 

A  4* 

A5* 

A  V 

0.00 

0.000 

1.406E-6 

1.812E-4 

2.212E-3 

6.913E-3 

7.446E-3 

2.989E-3 

0.25 

3.516E-7 

9.199E-5 

1.840E-3 

9.125E-3 

1.622E-2 

1.193E-2 

0.50 

2.335E-5 

1.012E-3 

8.684E-3 

2.534E-2 

3.1I3E-2 

0.75 

2.764E-4 

5.354E-3 

2.769E-2 

5.647E-2 

1.00 

1.615E-3 

1.920E-2 

7.004E-2 

1.25 

6.416E-3 

5.422E-2 

1.50 

1.997E-2 

Fig.  2.3.1:  Analytic  approach  to  computing  A  "exp. 

example:  Consider  the  second  table  in  Fig.  2.2.5.  The  value  0.01712  for  Aoexp  contains  no 
correct  digits.  Any  sixth  order  difference  of  a  polynomial  of  degree  five  is  zero 
(§2.7).  By  linearity 

A06exp  «  A<j fe , 

where  g  -  txp-ps  and  Ps  is  any  fifth  order  polynomial.  We  set 

the  first  terms  of  exp’s  Taylor  series.  For  this  choice. 

Using  g  instead  of  exp  we  get  Aoexp  *  Ao g  **  0.002989,  which  has  three  correct 
decimal  digits.  The  polynomial  Ps({)  which  dominates  the  information  in  the  left 
most  digits  of  exp({),  the  digits  lost  in  differencing,  is  removed  in  forming  #(£)• 
The  information  needed  to  give  Aoexp  accurately  is  retained  in  g({),  but  is  lost  in 
exp({)  because  too  few  digits  are  carried. 

There  are  many  cases  in  which  the  standard  divided  difference  scheme  (2.1.6)  works  very 
well  (Fig.  2.3.2).  Because  the  scheme  is  so  simple  and  computationally  fast  we  want  to  use  it, 
when  possible.  We  need,  then,  an  analysis  of  the  standard  formula  in  order  to  distinguish 
those  cases  where  we  may  wish  to  employ  it.  This  leads  to  criteria  for  deciding  when  to  use  it, 
rather  than  some  other  formula. 
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Divided  differences  of  expt0  by  the  recursive  formula 

A°expio 

A  'expio 

A2expio 

A5expio 

A4exp)0 

Asexp10 

A6exp10 

1.000 

4.472E+ 1 

1.000E+3 

1.492E+4 

1.668E+5 

!  .491 E  4-6 

1.113E+7 

ion 

1.218E+1 

5.449E+2 

1.219E+4 

1.817E+5 

2.031E+6 

1.818E+7 

ESI 

1.484E+2 

6.638E+3 

1.485E+5 

2.213E+6 

2.475E+7 

0.75 

1.808E+3 

8.089E+4 

1.808E+6 

2.696E+7 

1.00 

2.203E+4 

9.851E+5 

2.203E+7 

1.25 

2.683E+5 

1.200E+7 

1.50 

3.269E+6 

Fig.  2.3.2:  Recursive  scheme  on  exp|0(£)  =  «,l0{,  correct  Aoexp|0-  1.1 12E+7  to  4  digits. 
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2.4  An  error  growth  analysis  of  the  standard  divided  difference  formula. 

The  standard  divided  difference  algorithm  (2.1.6)  propagates,  and  may  magnify,  errors 
introduced  at  earlier  steps.1  Algorithms  which  exhibit  this  unfortunate  error  magnification  pro¬ 
perty  are  often  shunned  in  practice;  however,  (2.1.6)  is  just  too  attractive  from  the  point  of 
view  of  speed  and  simplicity  to  be  discarded  out  of  hand.  We  study  here  the  standard  scheme's 
error  behavior  and  obtain  error  growth  bounds.  This  analysis  provides  criteria  for  deciding 
when  to  employ  the  standard  scheme,  or  another  method,  to  compute  Ao/ 

We  analyse  the  error  propagation  in  a  typical  step  of  (2.1.6), 


Ad 7 


Ar'/'-Ad"1/ 

s„-Co 


(2.4.1) 


For  any  expression  g  let  Mg)  represent  its  computed,  or  "on  hand,"  value.  Employing  previ 
ously  computed  values  in  (2.4.1), 


M±of)  =  - 


/HA 


r '/•>  -/(Ad--1/) 


Define 


JHA0V>  =  A Sf  +  So'  ■ 


(2.4.2) 


8d'  is  the  absolute  error  in  expressing  Ad'/ by  ,/HAd'/>-  Then 


Ad'/  +  #u* 


(Ar’/^gr1)  -  (Ao'-1,/  +sd":) 
L-i  0 


and  so 


-  a  0y  + 


8n  —  1 C  rf  — l 

i  Qq 

in  —  40 


s 


It 

0 


sr'-sr1 

{.-Co 


(2.4.3  * 


(2.4.3)  represents  So  as  the  error  propagated  from  errors  in  the  /; — 1  -st  order  differences.  The 

tAlgorithms  which  magnify  previously  introduced  errors  from  step  to  step  are  often  re¬ 
ferred  to  as  "unstable."  This  term  is  commonly  applied  to  algorithms  for  the  numerical 
solution  of  differential  equations.  In  this  context  it  is  employed,  for  example,  in  texts 
by  Richtmeyer  and  Morton  (1967)  and  Gear  [1971], 
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growth  in  this  propagated  error  is  governed  by  (£„  -£o>"\  the  inverse  spread  in  the  abscissae. 
We  may  also  define 


fl(±oJ)  =  A0'7(l  +€0")  . 


(2.4.4) 


€o  is  the  relative  error  in  expressing  Ao'/by  its  computed  value  /(A0"/)-  Then 

a  r  '/•  ( i +« r 1 )  -  a,5'-  '/•  ( i +«o " 1 ) 


A  67  (1  +  «6')  - 


A67- ii  +  «r‘  + 


5«-5o 

Ar7(«r'-«r') 


A67«„-{o) 


and  so 


«1 


A6,~l/-(t|'~l-60"~l) 

A67-(C„-{0) 


(2.4.5) 


is  the  relative  error  in  /(A  6/)  propagated  from  relative  errors  in  n-l-st  order  differences. 

Expression  (2.4.5)  indicates  the  relative  error  may  grow  from  step  to  step  in  (2.1.6),  espe¬ 
cially  when  the  abscissae  £ o  and  £„  are  close.  This  relative  error  growth  is  equivalent  to  the  loss 
of  information  discussed  in  §2.2.  Such  growth  in  practice  may  nearly  approximate  the  upper 
bounds  on  error  growth  we  derive  in  §3.2. 

Insisting  on  small  relative  errors  is  often  inappropriate  in  divided  difference  computations. 
From  (2.4.5)  one  expects  a  large  increase  in  the  relative  error  when  |A6/|  is  small  compared 
with  |A6'_,/|.  However  a  large  relative  error  in  a  small  number  is  not  a  disaster  when  the 
absolute  error  is  small  relative  to  the  final  quantity  in  the  computation  in  which  the  divided 
difference  is  used.  Our  interest  here  is  accurate  computation  of  n-th  order  divided  differences; 
we  must  then,  at  least,  compare  the  absolute  error  with  an  appropriate  estimate  of  the  magni¬ 
tude  of  n-th  order  divided  differences  of  /.  Conclusions  regarding  the  bounding  of  the  errors 
expressed  in  (2.4.3)  and,  especially,  (2.4.5)  depend  on  the  particular  function  /  and  its  own 

divided  differences.  We  study  the  exponential  function  in  Chapter  3. 

/ 

example;  Large  relative  errors  in  small  numbers  are  not  always  disastrous.  Let  /-cosh. 

Using  four  digits  we  compute  the  following  differences.  The  entry 
A/cosh-5.970£-3  has  a  greater  relative  error  than  the  other  entries;  yet  subse¬ 
quent  table  entries  are  unaffected  by  this  error  because  the  number  5.970E-3  is 
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Table  using  4  digits 


0.005970 


Correct  values  to  4  digits 

4°cosh 

4  'cosh 

4  2cosh 

4  ’cosh 

- - 

3.762 

1.543 

1.555 

27.31 

-2.219 

0.005885 

8.613 

0.7392 

1.721 

0.1637 

Fig.  2.4.1:  Large  relative  errors  in  small  numbers  may  not  be  important. 


small.  In  four  digit  subtraction  0.005970- (-2.219)  forms  402cosh-0.7392:  the 
incorrect  rightmost  digits  play  no  role. 


2.5  The  divided  difference  as  a  function  of  its  abscissae. 

Our  approach  to  divided  differences  has  not  escaped  the  notion  that  their  formation  is  an 
operation  performed  on  entries  in  a  mathematical  table.  In  §2.3  we  illustrated  the  usefulness  of 
an  analytic  approach.  The  basis  of  this  approach  is  that  divided  differences  are  functions  of 
their  abscissae  and  may  be  treated  as  we  treat  other  mathematical  functions. 

To  aid  our  discussion  we  introduce  a  vector  notation  for  divided  differences. T  The 

sequence  of  abscissae  Z  —  |£o.  4i . £„l  is  conveniently  viewed  as  an  n+1  -tuple.  Hence  Z  is 

equivalent  to  a  vector  z  =  (£0, £i.  ■  ■  ■  ,  S„)  in  Cn+I  (or  Rll+I  for  real  abscissae).  We  speak, 
then,  of  a  divided  difference  function  A"/ being  defined  for  a  vector  z  in  the  same  sense  as  the 
function  /  being  "defined  on  the  sequence  of  abscissae  Z"  (§2.1).  Thus 

A*/(z)  =  U  "/)(£<>.£  . W-  (2.5.1) 

When  every  vector  z  in  a  region  of  C*+l  is  equivalent  to  a  sequence  of  abscissae  Z  on  which  / 
is  "defined,"  A  "/is  a  function  on  that  region.  Thus  when  defined, 

A"/  :  Cn+I  —  C, 

in  brief.  Our  new  notation  expresses  o-lh  divided  differences  of/,  A "f  as  functions  from  C"+l 
into  C.  The  value  of  this  function  at  the  point  z  €  Cn+I  is  A "/(:).  The  ordering  of  the  abscis¬ 
sae  is  suppressed  here. 

When  /  is  holomorphic  on  a  region  containing  the  abscissae,  A  "/is  holomorphic  in  each 
of  its  abscissae.  In  particular  for  each  /  *0. 1 . n. 

3  . .  ,  ..  a "/(z +  «',-£,)•*,)-  A"/(z) 

—-A'7(z)  -  im - - — - - 

Hi 

-  A"+l/(r,£,).  (2.5.2) 

The  vector  e,  =  (0.0 . 0,1.0 . 0)  is  the  /- th  coordinate  vector  in  C,H'1.  The  partial 

derivative  with  respect  to  £,  of  A  "/is  an  rr+I-st  order  divided  difference  with  the  abscissa  £, 
repeated;  this  is  indicated  by  (z. {,). 

fThe  first  notation  (A"/)({o.£i . {,,)  emphasizes  both  the  sequence  of  abscissae 

and  its  ordering.  The  second  notation  A <('/  merely  emphasizes  the  sequence  ordering;  it 
suppresses  reference  to  a  particular  sequence. 
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That  A"/  is  holomorphic  in  each  of  its  abscissae  separately  suggests  it  may  be  expanded  in 
a  series.  The  next  three  sections  develop  such  expansions. 

example:  By  differentiating  successive  divided  differences  we  can  show  that 

oo  (_r\k-n 

(A"exp)(0,£,£ . £)  -  . 

k-n  *  • 

Start  with 


ei  —  1  |  t  1  —  e  {  (— £);  1 


(A'exp)(0,£)  -  ,  -  e 4- — r —  -  e' 


A-l 


A:! 


With  £|-£, 


(A2exp)(0,£,£)  -  -—(A'exp)(0,£|)  -  -^r(A'exp)(0,£) 
o£  i  “4 


I- 1  '•  X— 2  K~ 


(Ar-l)!  ft!  ' 


*-2 


k\ 


Using  the  chain  rule  in  the  general  case. 


(A''+lexp)(0,  £ . £)  -  -~-(A''exp)(0.£|,£ . £)  -  —  ~ ( A  "exp) (0. £ . £) 

o£i  n  </£ 


*  ,r„  /!  *-Th  *! 


A  —  /i  —  I 


-  V  £ 


k ! 


(-4)*-"-1 
k\ 


This  may  be  compared  with  a  method  based  upon  the  standard  formula, 
(A"+lexp)(0,£,£ . £)  -  U"exP)({’{-  .  ~  (A"exp)(0.^.  -..{I 


j 


example:  Two-point  divided  differences.  Just  as  divided  differences  for  confluent  abscissae 
(which  may  be  referred  to  as  one-point  divided  differences)  reduce  to  a  special 
form,  divided  differences  about  two  repeated  abscissae  also  have  special  properties. 
Let  the  vector  r  =  .  £.  -()  consist  of  n  repetitions  of  the  two  data  points  £ 

and  -£.  Recalling  the  contour  integral  formula  for  divided  differences  (2.1.13),  for 
n  —0, 1,2,... 


_  liril  <«-{>*♦'(«  +  {)■  ' 


and 


1  r _ f((o)  do* _ 

2iriJc  (at-{)"+,(at  +  £)"+l 


The  first  In  abscissae  are  represented  by  z,  for  compactness.  For  each  n  define  the 
functions  b„  and  a„  by 


4.«)  =  +  a!7U-{)I  - 


o 


n+ 1 


and 


a  JO  s  A2"+l/(r,{,-{) 


.  :«+i 


/(at)  dot 


2niJc  (at-C)"+l(at  +  C)',+l 


These  functions  are  holomorphic  in  0 


2(«+l){-W{). 


and  similarly 

-fjjaM)  -2<*+l)£a,(+1(£). 

«£ 

The  functions,  then,  satisfy  the  recurrences 


b"(()  ~  2n('d(  6"'|(i) 


a"(S)  “  2n('  di  a"“|(£) 


and  can  be  defined  when  we  know  f>0(£)  and  a0(£)-  For  example  when  /— expr, 
expr(£)  =  eT(  with  t  >  0, 


b 0(£)  -  cosh(r£) 
a0(()  -  sinh(r£)/£ . 


Since 


a,M)  -A:"+1/(r,{,-{) 


A;'7(z.£)-A2"/(z,-£) 

2£ 


the  divided  differences  are  recoverable  from  &„(£)  and  a„(£).  The  values  b„(0  and 
a„(£)  yield  coefficients  of  a  Newton  expansion  of  /about  £  and  -£.  Note  that  when 
£-/ 1\  is  pure  imaginary,  both  &„(£)  and  a„(£)  are  real  for  any  /  such  that 
/(£)  —  /(£).  We  extend  this  example  in  §2.8. 
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2.6  The  divided  difference  table  as  a  function  of  a  matrix. 

The  entire  divided  difference  table  (§2.1)  of  /  for  the  sequence  of  abscissae 
Z  —  (Co. . . . {»)  may  be  expressed  as  an  (n+1)  x  (n+l)  upper  triangular  matrix' 


/(Co>  Ao/ 

X  If  ■ 

xu 

/<C.) 

A  ,7  • 

■  xr'f 

Ail)  ■ 

■  Aj •-*/ 

/(C„) 

Let  Zbe  the  special  («+l)  x  (/i-fl)  bidiagonal  matrix 

Co  1 
C.  1 
?2  1 

Z  -  (2.6.2) 

Cm— I  1 

Cm 

Opitz  [1964]  refers  to  Zas  a  "steigungsmatrix"  (ascent  matrix).  We  shall  call  Z  a  "step  matrix." 

The  same  conditions  on  /  imply  the  existence  of  both  the  divided  difference  table  Xf 
(§2.1)  and  the  Newton  polynomial  representation  of  /(Z)  (§1.2).  The  two  are  related  as  fol¬ 
lows. 

Theorem:  "The  divided  difference  table  is  a  matrix  function." 

A/-/(Z)  (2.6.3) 

proof:  The  Newton  polynomial  representation  of  /(Z)  is 

/(Z)  -  /(Co)  /  +  A«j/-(Z-C0/)  +  '  '  '  +  • 

*-o 

t  Az/is  written  when  the  sequence  of  abscissae  Z  must  be  emphasized.  Recall  that 
X  f,  no  superscript,  is  a  matrix  and  A '7  is  a  scalar  function. 
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Iff 

Because  Z-{*/,  each  k,  is  bidiagonal,  the  product  matrices  (Z  — <* /)  for  m  <  n- 1  have 

*-o 

II -I 

(0,/r)  element  zero,  while  the  (0 ,n)  element  of  fJ(Z  -{*/)  is  one.  Thus  the  (0,«)  element 

i- o 

of  /(Z)  is  A<f/  the  (0,/j)  element  of  A/  By  the  pattern  of  dependence  (§2.1),  the  choice  of 
0-th  and  n-th  abscissae  is  arbitrary.  Hence  equality  holds  between  every  element  of  ,/(Z>  and 
A./:  □ 

Parlett's  recurrence  (1.1. 5)  reduces  to  the  standard  divided  difference  scheme  (2.1.6) 
when  the  upper  triangular  matrix  T  is  replaced  by  the  step  matrix  Z.  This  provides  another  way 
to  establish  (2.6.3). 


Several  important  and  useful  consequences  follow  from  the  theorem. 


1. 


Function  of  a  Jordan  block.  When  the  sequence  of  abscissae  Z  —  lio.£o . Co)  >s 

confluent. 


/(Co)  /(Co) 

3?/’ (Co)  •  • 

7r/",(Co) 

/(Co) 

/(Co) 

/(Co) 

A  /- 


(2.6.4) 


I  /(Co)  I 

This  is  the  well-known  special  form  for  a  function  of  a  Jordan  block. 

2.  Multiplication  formula.  Let  the  function  fT  be  defined  by  fr(0  =/{rO ,  then 

A/r  -  /(tZ)  .  (2.6.5) 

3.  Scaling  abscissae.  Let  D  =  diag(l,T,T2 . r"),  a  diagonal  matrix,  and 

rZ  =  |r{0,  . . r{„),  then 

Az/r  ™  D'k,lf  D~X. 


proof:  Ai/r-/r(Z)-/(TZ)-/(DZrZr')-Z>  /(ZT)  £r'-£>  ArZ/  D-',  where 


(2.6.6) 


r(o  1 


^2  1 


T{„_,  1 


This  is  the  scaling  invariance  property  (2.1.10)  in  matrix  form. 

Special  functions.  Divided  difference  tables  of  certain  functions  inherit  some  appealing 
properties  from  the  functions  themselves.  For  example  when  /  —  f the  >th  power  func¬ 
tion  f '(()  =  ('  for  J-0. 1,2 . 


!}>+*-  Z'+A-  A  f'-Af*. 


Also  when  / -expT,  expr(£)  3  er{, 


AexpT+(r  —  p(r+,riz-  erZ  e"z  -  AexprAexp,r . 


(2.6.7) 


(2.6.8) 


tOur  divided  difference  notation  suppresses  variables,  so  clarity  demands  that  every 
function  have  a  name.  The  notation  f'  for  the  >th  power  function  is  used  by  Davis 
(19731. 


7 
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2.7  Divided  differences  of  polynomials. 

To  aid  in  the  development  of  series  expansions  for  A "f,  we  first  examine  divided 
differences  of  polynomials.  Let  pf  be  a  monic  polynomial  of  degree  j  given  in  factored  form* 

Pj(0  =  “  /»/-iU) •({-<*/_!)  (2.7.1) 

0 

for  any  ./- 0,1.2,...  {/>0(£)  =  l].  The  polynomial  p,  appears  in  the  >th  term  of  the  Newton 
expansion  (§1.2)  of  ./'about  the  sequence  A  s  {a0.aj.a2....), 

j-  0 

Pi  reduces  to  the  >th  power  function  ]J  when  all  the  a,  are  zero.  For  any  j  and  step  matrix  Z 
(any  sequence  of  n+1  abscissae),  the  matrix  function  theorem  of  the  previous  section  yields 

kPj-Pj(Z)  H  fl(Z-a,/).  (2.7.2) 

,-o 

example:  When  n  -  4  and  j  —  3, 

kPi  “  Pj(Z)  =  (Z -a0/)(Z-at|/)(Z-a2/) 

1 

£|-«,  1 

-  n  {j-®-  1 

.-0 

£j-«.  1 

S4~«, 


in  non-factored  form,  the  linearity  property  (2.1.8)  yields 

/•0 


tFor  a  polynomial  p  - 

-AT- 
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(4o  —  ao)(£o~a|)(£o  —  ®2^  ££(£*  ~ ~«») 

A  — Oi— 0 

(£t  “ao)(£i  —  <*i)(£i  —  <*2) 


£(£.-«.) 

t~0 

£  £(£t+i  “®t+iKC,+i  —  a,) 

k-0i-0 

({2-ao)(C2-«i)({2~«2) 


1 

£({,+i-at,) 

i-0 

lX({*+2“aH|)({i+2”ai) 
*  -0/-0 

(£3~«oH{3-“l)(Cj-«2) 


0 

1 


I  k 


X£^*+i-a*+l^i'+3~a^ 

i-0 1-0 


(£4  ~  «o)  (£4  ~  « |)  (£4  —  «2^ 


In  particular,  the  0-th  (top)  row  of  lp}  is 
A0P3 "  (£o“«o)(£o~«i)(£o~“2) 


&0P3  “  (£o~«o)(£o~«i)  +  (£o~«o) (£1  “<*2)  +  (£1 —  “i)(£i  —<*2) 


^oP3  "*  (£o~<*o)  +  (£i  —  ®i)  +  (£2  —  <*2^ 

A<jV.3  “  1 

Aq pj  ™  0 


For  general  n  and  j  the  (0,n)  element  of  ip,,  that  is  A <?/>/,  is  the  (O./r)  element  of  the 
matrix  product  ( Z-a0/)(Z  ~a\ /)  •  •  •  (Z  -o,_|/).  The  following  formulas  can  be  verified  by 

k  k 

actually  writing  out  the  products.  We  freely  use  the  convention  that  JJs,  =  1  and  £s(  s  0 

<•/  /»/ 

when  k  <  j,  where  s,  represents  some  expression. 


y 
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In  the  special  case  of  the  >th  power  function  f7,  expression  (2.7.4)  reduces  to 

AoT-  I  M' 

A0+*l+  •  ■  ■ 

k,Z  0 

This  is  a  well-known  symmetric  polynomial  formula  for  the  divided  difference  of  power  func¬ 
tions  [Milne-Thomson,  19331.  When  n  —  1  the  first  divided  difference  of  pj  obtained  from 
(2.7.4)  is 

A:*0  /• 0  i«-0 

-  glfltto-.,)  fl  •  (2.7.5) 

k~Q  /■» 0  (•Ir  +  I 

A  simple  recurrence  for  computing  divided  differences  of  the  functions  p,+ 1, 
j  “  “l . 0. 1 . 2. ....  may  be  developed  from  expression  (2.7.2).  We  begin  by  writing 
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Lpl+l-  kpj  (Z-ajl).  (2.7.6) 

By  writing  out  the  right-hand  matrix  product,  the  (O.k)  element  of  A pi+ ,  is 

AoP,+i  -  ({*  -a,)-A0*P/  +  (2.7.7) 

for  any  *  -0. 1 . n.  All  the  elements  Afipj+t  for  0  ^  k  ^  n  comprise  the  0-th  (i.e.  top)  row 

of  the  matrix  Api+i. 

Formula  (2.7.7)  is  a  recurrence  in  k  and  j.  To  see  this,  we  replace  the  index  j  in  (2.7.7) 
by  j+k  where  0  <  k  <  n,  and  still  j  -  - 1 , 0. 1 , 2 .  That  is, 

AoPj+k+\  ™  ({*  —  a)+^)-Aop;+t  +  Ao  'p./**  (2.7.8) 

is  the  (0,*)  element  of  the  matrix  Apy+*+|.  Thus  for  fixed  >,  varying  k  in  (2.7.8)  has  us  look¬ 
ing  at  elements  from  the  top  row  of  different  matrices. 

example:  Let  j  —  2,  then  for 

k-0,  Ao*P/+a+i  ”  AoPj  is  the  (0,0)  element  of  A  p3; 

k  - 1.  AoV;+*+i  “  A0V4  is  the  (0, 1)  element  of  A  p4; 

k  -  2,  AoV,+*+i  “  A0P5  is  the  (0, 2)  element  of  Aps; 

and  finally  for 

k  -  n,  AoP/>a+i  ”  Ao'p„+}  is  the  ( O.n )  element  of  Ap„+i . 

Since  AoP*-l  for  ail  k,  the  elements  AoP/+*+i  are  known  for  /  -  - 1 ;  so  all  the  AoPa+i 
are  defined  by  (2.7.8)  [we  define  A-'p*  =0  for  any  k).  Thus  all  the  AoP/+t+i  are  computable 
for  j  —  0,  and  recursively  for  any  J  >  0  as  well.  This  procedure  is  summarized  in  Algorithm  1 , 
and  its  first  few  steps  are  illustrated  in  Fig.  2.7.1.  Note  that  if  we  want  all  the  top  row  elements 
of  the  table  Apm+h  m  ^  rt,  one  element  appears  in  each  step  of  the  algorithm  from  j  -  m-n  to 
j  -m.  AoPm+i  appears  first  (step  j  -  m—n)  and  AoPm+i  appears  last  (step  j  —  m).  Each  >step 
of  the  algorithm  requires  n+\  multiplications.  Three  storage  n + 1  -vectors  are  needed:  one  to 
hold  the  abscissae  one  to  hold  the  >th  level  results  AoPju-,+i  (the  results  for  level  j— l  may 
be  overwritten),  and  one  to  hold  the  n+I  currently  active  <*,,  namely  c*,,  a)+| . a;+„. 
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Algorithm  1:  Recursive  computation  of  A &Pt+j. 

1.  Initialize  A  oft  - 1  for  each  k  —  0, 1. . .  .  ,n 

2.  For  7 -0.1.2.  ... 

Ao/>/+*+i  ”  (£*  — «*+>)' A^ft+>  +  Ao_,ft+y*  "* 0,1... 

.  .  n. 

(A  kf  s  0  for  k  <  0,  any  function  ,/) 

Initialize 

A oPo  —  Aoft  —  4<Jp2  =  I 
For  7  -  0, 

A$p i  "  (£o— ao)‘AoPo  +  Ao'po “  £o~ ®o 

A<Jp2  *“  (£i  ~ a()-Ao^i  +  AoPi  “  (£i  ~ «|)  +  (£o~ “o) 

A02/>j  *■  (£:  — a2)'^(?P2  +  A0P2  “  (£2“ “2)  +  (£1 —  «i)  +  (£o~ «o) 

For  7-1, 

A0P2  “  (£o~®i)  b§P\  +  A<j"  ]P\  “  (£o— «i)(£o~«o) 

A0P3  ”  (£1  ~«2)'AoP2  +  &8pi 

—  (£|  —  aj)(£i -ai)  +  (£1 —  «2^£o— ®o)  +  (£o  —  ai)(£o—ao) 

A^ft  ”  ((2~a3)'&rfpi  +  A^J 

—  (ii  —  aJiCi  —  ati)  +  (£2" <*?)(£i  ~ai)  +  (£:~ °3^£o  —  0io) 

+  (£1  ~<»2H£i  ~0|)  •+■  (£1 —  a2^£o~tto^  +  (£o_«i)(£o“ao) 

Fig.  2.7.1  :  First  couple  of  steps  of  Algorithm  1  for  n  -  2. 

A  companion  algorithm  for  computing  the  n-th  column  of  Apf+i  also  exists.  Obtaining  it 
merely  requires  rewriting  (2.7.6)  as 


§2. 


/-i 

A/»;+ 1  ••  (Z  —  a,/)-JJ(Z  —  a,/)  *  (Z  —a,l)  kpj  (2.7.9) 

,-a 

and  following  the  same  approach  as  before.  Again  just  one  element  of  the  n-th  column  of  a 
particular  matrix  is  computed  at  each  step.  The  first  few  steps  of  the  algorithm  are  illustrated  in 
Fig.  2.7.2. 

Algorithm  2:  Recursive  computation  of  A  *_*/>*•»■/• 

1.  Initialize  A "  1  for  each  k  -0. 1 . n 

2.  For  jm  0,1,2, ... 

A*_aa+/+i  ■  -<*k+,)-b«-kPk+,  +  ^n-Ui Pk+,-  *"0,1 . n. 
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Initialize 

^  iPo  “  *iP\  ”  S  I 

For  j  -  0, 

^iP\  “  (i2~ao)^2PO  +  "  £2“ «0 

^i'P2  “  (£1  —  c*  1 ) •  A i1  />  1  +  <l2/>i  ”  (£1  —  at  1 )  +  (£2  —  a0) 

&$Py  “  (Co~at2)'A<}p2  +  ^i'P2  ”  (€0  — «2>  +  (£i  — ®|)  +  (£2“ ®n) 

For  j  —  1, 

^2/,2  ”  (£2  —  ®l)‘^2^l  +  Af'*l  “  (£2  — tt))(C2“  ao) 

&\Pi  ”  (Cl  —  <*2) l*/?2  +  ^2°P2 

“  (£l  — <*2)(£l  — al)  +  (C|-«2)(£2~«o)  +  (C2~al)(£2~ao) 

^$P*  ”  (^o  —  «3>*A(5 |/7j  4-  A)1/?.! 

*■  (£0** <*.i)(£o"“a2)  +  (£o“aj)(£i“«i)  +  (£o  — a3)(C2  — ao) 

+  (£1 -a2>(£i  — a|)  +  (£  1  —  at2)  (£2  —  “o)  +  (£2-ai)(C2~«o) 

Fig.  2.7.2:  First  couple  of  steps  of  Algorithm  2  for  n  —  2. 


/ 


/ 
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2.8  Series  expansions  of  A"f. 

The  matrix  function  theorem  for  divided  difference  tables  leads  directly  to  series  represen¬ 
tations  for  divided  differences.  For  example  suppose  that  on  the  disk  |p>  |£ -a|),  / 

has  a  Newton  expansion  about  the  sequence  of  expansion  points  A«  {a0.at.a2....).  That  is 


(2.8.1) 

A-0 

*-l 

on  D,„  where  />*(£)  —  JJ  (£—«*,).  Appendix  A  sufficient  conditions  are  presented  for  the 

J-0 

existence  of  such  an  expansion.  Under  the  same  conditions  the  matrix  function  f(A )  has  a 
Newton  expansion  when  ail  the  eigenvalues  of  A  lie  in  D/r  Thus  when  the  data  points 
Z-{£0.£i . £„)  lie  in  D,„  the  divided  difference  table  has  the  Newton  expansion 

A/-/(Z)-  £  A* J-pk(Z):  (2.8.2) 

*-0 


In  the  vector  notation  introduced  in  §2.5,  Z  C  D;,  is  equivalent  to  a  vector 
-”(£o.£i.  •  •  •  .£«)  in  D"+l  CC"+I.  The  previous  section  leads  us  to  examine  the  (0 ,n)  ele¬ 
ment  of  the  table.  The  result  is  summarized  in  the  following  theorem. 

Theorem:  Newton  expansion  of  the  divided  difference  function.  Suppose  /  has  a  Newton 
series  on  D,,.  Then 

A"/-lAay-A>*  (2.8.3) 

*-1 

over  all  D"+\  where  pk(0  -  n  (£-«;). 

./•0 


The  important  point  is  the  identification  by  (2.7.2)  of  the  (0 ,n)  element  of  pk(Z)  with  Ao 'pk. 

A  Taylor  series  expansion  formula  for  A"/  is  an  immediate  corollary.  Recall  that 
A«(/“/<*)(a )/k\  in  the  confluent  case. 

tThe  reader  is  asked  to  distinguish  between  the  divided  differences  A ^  forming  the 
series  coefficients  which  have  abscissae  in  A,  and  the  elements  A"/  of  the  divided 
difference  table  which  have  abscissae  in  Z.  Brent  [1973)  presents  a  simple  Taylor  ex¬ 
pansion  for  the  divided  difference. 

\ 


J 
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Corollary:  Taylor  expansion  of  the  divided  difference  function.  Suppose  /  is  holomorphic 
on  a  region  D  containing  the  point  a.  Let  D,,-{£|p>  |£-a|)  be  the  largest  open  disk 
such  that  D,,  C  D.  Then 

IV  -  I^"a)  l"TaA  12.8.4) 

k -n 

over  all  D*+l,  where  ti(£)  =  t‘(£ -a)  -  (£  ~a)k. 


proof:  Because  /  is  holomorphic  on  D,„  it  has  a  Taylor  expansion  about  a  for  all  £  €  D„.  The 
theorems  in  Appendix  A  establish  that 

f(Z)  -  I  L-:~}k(Z-aI) .  □ 

k-0  *  ■ 

Formulas  (2.8.3)  and  (2.8.4)  suggest  ways  to  compute  divided  differences  for  perturbed 
abscissae  when  the  unperturbed  divided  differences  are  available.  The  computation  of  divided 
differences  by  (2.8.4)  for  functions  such  as  exp,  sin,  and  cosh  is  quite  straightforward  since  the 
Taylor  coefficients  are  easily  obtained.  Functions  such  as  log  and  may  also  be  treated;  how¬ 
ever,  care  is  required  to  ensure  that  we  use  a  series  representation  whose  circle  of  convergence 
contains  all  the  data  points. 

The  algorithms  of  §2.7  in  combination  with  (2.8.3)  lead  to  a  method  for  computing  A0‘/, 
0 <  k  <  n,  when  we  already  know  the  coefficients  Ao^s/3,,  /«■  0.1,2,...,  of  / s  Newton 

in 

expansion.  Let  s„,  =  £)3/Pf  be  the  partial  sums  of  the  Newton  expansion  (2.8.1)  of  /  so 

i-o 

s,„  —  /as  m  —  oo.  Then  by  linearity 

m 

Ao ‘sm  * 

l-k 

and  by  (2.8.3)  — 1  A  o' /  as  m  —  °°  for  any  k.  The  following  algorithm  computes  Jko/for  all 

k -0,1,  .  .  .  ,n  by  forming  the  partial  sums  A 0*Sw  for  m-j+k+l.  One  additional  term  is 
added  to  Ao*/+*+i»  for  each  k ,  at  each  /step. 
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Algorithm  1 :  Newton  expansion  of  Ao /. 

1.  Initialize  A0ap*  ”  1 .  Aos*  *ft .  for  each  k  »0. 1 . n 

2.  For  j -0.1.2. ... 

AoP)+a+!  ■*  (£*  ~ a/+A)'AoP/+*  +  Aoa-iP/+a 

Ao*5/+a+i  “  A0*j/+a  +  ft+*+f±<fP/+,i+i .  for  each  k- 0.1 . n. 

Exclusive  of  the  coefficient  evaluations,  the  scheme  requires  2w+2  multiplications  per  >step. 
Initialize 

AoPo  “  A(jpo  ”  AoPo  =  1 
Ao-o”/3oi  AoS|  —  ft,  A0S2  “  ft 

For  y  -  0, 


A0°Pi 

-  (Co- 

■ao)‘AoPo 

+  AiT'po  * 

■  ?o-«o 

Ao°/  : 

Ao°^i 

"  Ao°Jo  + 

ft'A<?,Pi  ** 

ft  +  0l((o~ao) 

A0P2 

-  (i,  - 

■a|)AoPi 

+  A^Pi  — 

((i—aj)  +  ({0-«o) 

An'y  ■ 

Aoi2 

-  Ao'rt  + 

ft'A0V:  " 

ft  +  02'Kil  ~  «i)  + 

^0~ao)l 

AoVt 

-  Uz - 

■a2)‘A(Jp2 

+  A0P2 " 

(<2-02)  +  (i|-«i) 

+  (Co~«o) 

Atf/ 

AoJj 

-  A02J2  + 

^•A^j  - 

ft  +  + 

(£i-«i)  + 

Fig.  2.8.1:  First  step  of  Algorithm  1  when  /r-2. 


There  is  also  a  companion  algorithm  which  computes  the  w-th  column  of  the  matrix  A/ 


PT* ' 
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Algorithm  2:  Newton  expansion  of  A  *_*/. 

1.  Initialize  -  1 .  A  “0k .  for  each  k  —  0, 1,  .  . . ,  n 

2.  For  7-0, 1.2.... 

kh-kPj+k  + 1  ”  ({«-*  —  <*j+k)'&!t-tPi+k  +  ^h-k  +  lPi+k 

A,t-*s/+*+l  -  kZ-kSf+k  +  0,+k*\'k^-kPi+k+\ .  for  each  *-0,1 . n. 


example:  Two-point  divided  differences  (cont.  from  §2.5).  The  Newton  expansion  of /about 
a  sequence  of  two  repeated  points  A  —  [a.  —  a.  a.  —  a, ...)  is 

/({)  -  /(a)  +  (A'/Ka.-aMS-a)  +  (k2f)(a,-a,a)((2 -a2) 

+  (k3/)(a,—a,a,—a)((2  —  a2)((—a)  +  ■  ■  ■  . 

We  may  also  expand  about  the  rearranged  sequence  -A-{—  a,cr,-a,a, ...), 

/(£)  —  /(— «)  +  (A '/)(-■< a.aMC  +  a)  +  (k2/)(— a, a,—  «M{2  —  a2) 

+  (k}/)(—a,a.—a.a)-((2-a2)((+a)  +  •  •  • 

Recalling  the  definitions  of  the  functions 

6„(a)  -  y{(A2"/)(a, -a,  .  .  .  ,a)  +  (A2"/)(— a, a . -a)) 

and 

a„(a)  —  (A2"+l/)(a,-«,  .  .  .  ,a,-a) 

from  §2.5,  the  average  of  the  two  expansions  is 

/(£)  -  bo(a)  +  +  6|(a)-(£2-a2)  +  a  |(a)(£2-a2)£  +  •  •  . 

We  remarked  earlier  that  when  a  —  irt  is  pure  imaginary,  a„(a)  and  b„{a)  are  real 
for  any  /such  that  /({)  - /(£).  Hence  the  expansion 

/({)  -  M/tj)  +  PotivH  +  b\(iriMl2  +  T)  +  a|(n?)  ({J  +  V)£  +  '  '  ' 


is  entirely  real  when  (  is  real. 
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1 


example:  Two-point  divided  differences  of  the  exponential.  When  /-expr,  expT({)  =  i 


with  r  >  0,  the  coefficients  bja)  and  a„( a)  satisfy  the  recurrences 
ra(,_,(a) 


b„(a )  - 


2/7 


a„(a)  - 


rb„-\(a)  -  (2/7-l)aB_l(a) 


2  not1 

To  show  this,  recall  that  (§2.5) 


b0(a)  “  cosh(ra) 


a0(a)  «•  sinh(T<*)/a 


and 


.  /  >  1  d  ,  ,  ,  rsinh(ra)  ra0(a) 

- £ - 2~ 


«,(«)  -  -L.-f  ao<„)  .  IS5a?(~I-ilnh(r.)/« 

2a  da  2a2 


rb0(a)  —  a0(a) 


2  a2 


In  general  (suppressing  a  for  compactness), 


1  (  -  (2/>-l )a„_i). 


"  2/ta2  2(n—\) 

Differentiating  we  obtain 


-^l^-^-DaV.) 


"  2na'  2(/7— 1) 


a 


or 


2(/7+l)a-a,l+|  -  —  -  T{raia^-l-2/7(2/i— l)a-a„)  - 

2na  a 


“  —  ( ■  -  (2/7-1  )a„ ) , 

a  2  n 


where  we  have  used  that  a'w_t  -2/ra-a„.  Thus 


„  1  !  Tfl»-i  ,  r&„- (2/7+1  )a„ 

2(fl+l)a2  2*  2(*+l)a2 


This,  along  with  a  similar  result  for  b„+ 1,  establishes  the  recurrences. 

For  example  when  a  -  5/  and  r  -  1 ,  we  obtain1 

bQ(5i)  -cos(5)  -  2.837£— 1 
a0(5i)  -sin(5)/5  -  -1.918£-1 
*,(5/)  -  -9.589£-2 
<7,(5/)  -  -9.509£-3 
b2(5i)  -  -2.377£-3 
<j2(5/)  -  6.737£-4 
b2(5i)  -  1 . 1 23 £ — 4 
<7j(5/)  -  3.830£-5 
b<(Si)  -  4.788£— 6 
<74(5/)  -  7.792£— 7. 

When  f- 

expd.Sir/)  -0.2837  -  0.1918  x  1.5ir/  -  0.09589  x  1(1.5#/)*- (5/)2]  +  •  •  • 

-  -1.000/ ,  as  it  should. 

example:  Suppose  we  want  to  compute  A02exp  for  the  slightly  perturbed  abscissae  Co-5.01/. 
C i  -  -5.01  /  and  Cj "  4.99/.  Let  us  follow  the  steps  of  Algorithm  1 . 

Initialize 

A0°l  -  A«j(t-5i)  -  <M(£  +  5/)  -  Ao2(£2  +  25)  -  1 
A0°Jo  -  0.2837 .  - -0.1918.  \$s2  -  -0.09589 

For  y-0, 

A0°(C-5/)  —  (5.01/  —  5/)  x  1  -0.01/ 

A0°(C  +  5 z)  -  (5.01/  +  5/)  x  1  -  10.01/ 

A0°£  -  5.01/ 

A0°s,  -  0.2837  -  0.1918  x  5.01/  -  0.2837  -  0.9608/ 

fMost  of  the  numerical  examples  here  were  done  on  a  pocket  calculator  which,  unless  a 
particular  working  precision  is  specified,  carried  more  digits  than  are  shown. 


I 
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Ad(£2  +  25)  -  (—5.01/ +  5/)  x  1  +0.01/  -0 
Ao:s2  -  -0.1918 

A0j({2  +  25)({-5/)  -  (4.99/  —  5/)  x  1  -+  0  —  —0.01/ 

A02(£2  +  25)(£  +  5/)  -  (4.99/ +  5/)  x  1  +  0  -  9.99/ 

Ad({2  +  25)£  —  4.99/ 

A,fo  -  -0.09589  -  0.009509  x  4.99/  -  -0.09589  -  0.04745/ 

For  7-1, 

A0°(?2  +  25)  -  (5.01/ +  50x0.01/  -  -0.1001 

A0°s2  -  0.2837  -  0.9608/  -  0.09589  x  (-0.1001)  -  0.2933  -  0.9608/ 

A<J(£2  +  25)(£-5/)  -  (-5.01/ -50  x  0  -  0.1001  - -0.1001 
Ao (£2  +  25)(£  +  5/)  -  (-5.01/  +  5OX0  -  0.1001  -  -0.1001 
Aj(£2  +  25)4:  -  -0.1001 

Ao-Jj  —  —0.1918  +  0.009509  x  0.1001  - -0.1908 

A(?({2  +  25)2  —  (4.99/  +  5/)  x  (—0.01/)  -  0.1001  -  -0.0002 

Ao2J4  -  -0.09589  -  0.04745/  -  0.002377  x  (-0.0002)  -  -0.09589  -  0.04745/ 

The  algorithm  may  be  continued  for  y-2,3 .  To  four  figures,  the  correct  values 

are 

A0°exp  -  0.2932  -  0.9560/ 

A(jexp  -  -0.1908 

A«?exp  -  -0.09589  -  0.04745/ 

Note  that  the  conjugate  pair  divided  difference  (A 'exp) (5.01/, -5.01/)  is  real. 
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2.9  Computing  divided  difference  tables. 

We  have  presented  essentially  two  very  different  methods  for  computing  divided 
differences.  The  first  was  the  standard  algorithm  (2.1.6)  which  is  simple  and  fast,  but  may 
magnify  errors  from  step  to  step.  The  second  was  the  more  complicated  series  algorithm  of 
§2.8. 

Because  Taylor  series  coefficients  are  most  easily  obtained,  the  series  algorithm  is  most 
easily  applied  for  a  single  expansion  point.  When  the  abscissae  are  closely  enough  clustered 
about  this  expansion  point,  the  series  is  rapidly  convergent.  Hence  the  series  algorithm  need 
be  computed  for  only  a  few  steps  to  obtain  divided  differences  with  small  error.  This  is  pre¬ 
cisely  opposite  to  the  case  for  the  standard  algorithm,  where  in  §2.4  the  error  magnification  was 
seen  to  depend  inversely  on  the  separation  of  the  abscissae. 

A  general  purpose  algorithm  for  computing  divided  difference  tables,  then,  will  be  a 
hybrid.  Each  algorithm  above  will  be  used  where  it  is  best  suited,  with  primary  consideration 
given  to  speed  and  accuracy  of  computation.  The  question  is  then  to  decide  which  method  to 
use  for  a  particular  element  of  the  table.  This  is  the  prime  topic  of  Chapter  3  where  divided 
differences  of  the  exponential  function  are  discussed. 

The  series  algorithm  presented  in  §2.8  computes  only  one  row  (or  column)  of  the  divided 
difference  table.  It  could  have  been  written  in  matrix  form  in  order  to  give  the  entire  table  at 
once.  This  is  equivalent  to  applying  the  given  algorithm  on  each  row  of  the  table. 

Such  repealed  applications  of  the  series  algorithm  is  not  necessary,  however.  After  one 
application  of  the  algorithm,  the  0-th  row  of  the.  divided  difference  table  is  obtained.  We  now 
have  sufficient  information  to  fill  out  the  remainder  of  the  table,  row  by  row.  by  running  the 
standard  scheme  (2.1.6)  backwards. 

Backfilling  the  divided  difference  table.  When  divided  differences  A0‘/  for  k  -0. 1 . n 

are  known,  the  remainder  of  the  table  may  be  obtained  by  computing  successively  for 
/  “1.2 . n 

A,*/  -  (U*-{,-i)-A,*-V/+  A*-,/.  (2.9.1) 


k  “0. 1 . n— i. 


Fig.  2.9.1:  Possible  order  of  backfilling  using  (2.9.1). 


As  in  §2.4,  we  check  how  errors  are  propagated  in  one  step  of  the  backfill  algorithm,  say 
to  compute  A A<T'/+  (£„-£<Mo/-  The  absolute  error  is 

Sf'  -  So"1  +  (£»-Co)«o;  (292) 


the  relative  error  is  given  by 


ir1/ 


Ad7 

kr'f 


a  „-co)«o'. 


(2.9.3) 


The  absolute  error  growth  is  governed  by  |{„  -£ol-  When  Ad'-1/  and  A \f  are  of  comparable 
magnitude,  the  coefficient  of  «d' .governs  the  growth  of  the  relative  error.  In  both  (2.9.2)  and 
(2.9.3)  the  expressions  governing  error  growth  are  essentially  inverses  of  those  governing  error 
growth  in  (2.4.3)  and  (2.4.5),  respectively.  Thus  the  backfill  algorithm  is  most  attractive  pre¬ 
cisely  when  the  series  algorithm  is  most  attractive  and  the  standard  scheme  is  not. 
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3.  Divided  Differences  of  Exponentials 


3.1  Special  formulas  for  exponential  divided  differences. 

The  ideas  presented  in  the  last  chapter  are  illustrated  by  considering  the  exponential  func¬ 
tion.  Because  the  exponential  function  is  entire,  all  formulas  and  algorithms  discussed  so  far 
are  applicable  for  any  abscissae.  Further,  the  special  properties  of  the  function  permit  useful 
simplification  of  our  previous  formulas.  In  addition,  results  obtained  for  the  exponential  may 
be  modified  to  cover  related  functions,  such  as  sin  or  cosh,  by  means  of  the  linearity  property 
(2.1.8)  of  divided  differences. 

The  behavior  of  exponential  divided  differences  under  a  constant  shift  in  the  abscissae 
illustrates  a  useful  simplification  of  the  translation  invariance  property  (2.1.9).  It  is  convenient 
to  consider  the  more  general  function  expT  with  scaling  parameter  r,  that  is  expr(£)  =  eH. 

Translation  property  of  exponential  divided  differences.  Let  z  be  a  vector  whose  elements 
are  data  points  (§2.5).  Then  for  any  constant  a, 

A"expr(z+av)  -  era  l"expTh)  (3.1.1) 

where  the  constant  vector  w  —  (1,1 . 1). 


j 


1 


/ 


It  is  clear  from  (3.1.1)  that  no  generality  is  lost  when  we  restrict  attention  to  sets  of 
abscissae  with,  say,  £o*0  or  C„““Co-  In  the  latter  case  the  first  divided  difference  simplifies  to 


A<jexpT 


e"T(o_er<o 

-Co- (Co) 


sinh(T{0) 

To 


In  general  for  non-centered  abscissae  any  first  divided  difference  of  the  exponential  can  be  writ¬ 
ten  as 


A(JexpT 


sinh(ri|/) 

* 


(3.1.2) 


where  to  =  ({,  +  Co)/2  and  <1/  s  (C  i  -  Co)/2. 

The  integral  representation  formula  for  divided  differences  (2.1.12)  acquires  a  simpler 
form  when  the  parameter  r  is  non-negative.  We  have 
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I  TJ  T*-l 

<M'exp,  -  ff  ■  ■  ■  f  -j—txptllo+  (Ci_Co)t|+  •  •  •+  (C„~{«_i)r„J  t/r„  ••  ■  rfr, 

0  0  o  “4 

l  ri  r»-i 

-  ff  •  •  J  T"exp{r{0+({|-{o>TT|+  •  •  ■  +({„-C„-|)rrMJrfT„  •  •  •  rfr, 

0  0  0 

by  the  definition  of  expr.  The  change  of  variables  ■*  rr ,  for  j  -  1 , 2 . n  yields  the  alter¬ 

native  expression 

T  "l 

^d'expT  ”  J*J*  /  explT{0+({,-io)(T,+  •  •  • +(£„-£B_|)<r„]</<7„  •  •  •  t/cr,  .(3.1.3) 
0  0  0 

We  recognize  that  this  is  a  recurrence  for  Ad'cxp,,  namely 

r 

Ao'expT  -  eiof  i"_,expir  da 

o 

where  cr  =  <r  i.  By  the  symmetry  property  (2.1.3),  the  ordering  of  the  abscissae  is  arbitrary;  we 
may  replace  Co  by  any  0  <  /  <  n.  To  deal  with  such  cases  we  define 

l/ir'expr  =  U-'-'exPrHCo.C . . £„).  (3.1.4) 

the  tf-l-st  divided  difference  with  the  Ath  abscissa  omitted.  (3.1.5)  summarizes  the  formula. 

Recursive  integral  formula  for  l^exp, .  For  any  r  ^  0  and  index  /  —  0, 1 . /?, 

r 

Ad'expT  -  eTi'f  e'^'  XC'j'sxPvdcr  .  (3.1.5) 

o 

This  result  will  prove  useful  in  the  next  section.  In  addition,  one  recognizes  that  formula 
(3.1.5)  is  a  convolution, 

(exp{  *  l/;7'exp)(r) , 

where  A  ",7 ‘exp  is  treated  as  a  function  in  r.  The  correspondence  is  obvious  from  the  convolu¬ 
tion  formula,  with  gfal-A/'T'exp,,, 

r 

(f*g)(r)  -  f  /{t  —<r)  g{<r)  d<r  . 

o 
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3.2  Real  exponential  divided  differences. 

Exponential  divided  differences  for  real  abscissae  are  positive  and  increasing  functions  of 
their  abscissae.  These  properties  permit  derivation  of  both  upper  and  tower  bounds  which  show 
how  the  divided  differences  and  the  error  growth  in  the  standard  formula  depend  on  (/)  the 
spread  in  the  abscissae.  («)  the  order  n  of  the  difference,  and  (///)  the  parameter  r. 

in  this  and  the  next  four  sections  we  consider  exclusively  divided  differences  of  the  func¬ 
tion  /— expr,  with  parameter  r  >  0,  for  real  sequences  of  abscissae  X  =  {£0.£i . £,,!•  All 

such  divided  differences  have  two  properties  which  characterize  them. 

Theorem  I:  For  all  r  >  0  and  n  >  0,  A"expT  is 
(/)  positive, 

(ii)  strictly  increasing  in  each  abscissa  i  -0, 1 . n. 


t 

i 

l 


!  ’ 


proof:  The  result  is  almost  immediate  from  the  recursive  integral  formula  (3.1.5), 

T 

A<j'expr  -  eT(  J e  ^'-A/ir’expff  da . 
o 

All  0-th  order  real  exponential  divided  differences  are  positive  (A°expr(f)  -eT*].  The 
recurrence  implies  all  first  order  differences  are  also  positive,  and  hence  by  induction  all  /Mh 
order  differences  are  positive  for  any  n.  For  (/>), 

_  r 

-r--Aoexpr  ™  [(T  —  a)eir  "i('-\",Jlsxp^da  >  0, 

O?/  1) 

since  the  integrand  is  positive.  □ 

The  recursive  integral  formula  also  provides  an  easy  way  to  develop  expressions  relating 
divided  differences  of  orders  n— l  and  n. 


•  " 

/ 
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Theorem  2:  Suppose  0  <  ( ,  4  y  for  every  abscissa  0  ^  i  ^  n.  Then  for  each  /  there  ex¬ 
ists  a  Z  €  1/3.  y]  such  that 


A/;,expr  -  +  —  )l0"exPr- 


(3.2.1) 


proof:  By  the  translation  (3.1.1)  and  scaling  invariance  (2.1.10)  properties, 

r 

l"exp(r(ar-fi/)l  -  r-'V-Tf  Ad'expT  -  r~“eT<f,~(’f  e~'r('-A<"j'ex p,r  c/<r 

it 

for  any  i  -0. 1 .....  /r  and  Z-  Differentiating  with  respect  to  r  yields 

~A"cxp[r (*-£«)]  -  T-f“T*|(f,-f--J-)-lo"hx pT  +  A/'j'exp,) .  (3.2.2) 

Every  element  of  the  vector  x-@uis  non-negative,  and  so  l"exp{r(x-/3u)l  is  increasing  in  r. 
Similarly,  every  element  of  x  —  yu  is  non-positive  and  A"exp[r(jr -yu)]  is  decreasing  in  r. 


Hence 


-~A"exp[r(x~0u) )  >  0  >  -f-A"explr(x~yu)]  ■ 
a  t  at 


so  for  some  Z  6  [/3.  y],  the  derivative  is  zero.  □ 


A  plethora  of  upper  and  lower  bounds  on  Ay'exp,.  can  be  derived  from  the  simple  expres¬ 
sion  (3.2.1)  by  choosing  particular  values  of  Z  and  i.  The  two  simplest  follow  by  choosing  Z  as 
one  or  the  other  of  the  end  points  of  the  smallest  interval  containing  the  abscissae.  Note  that 
equality  holds  when  the  abscissae  are  confluent. 


Corollary  1:  Lower  bound  on  A<j'expT.  If  ^  for  each  /  «-0, 1 . /r,  then 


Ad'expT  >  — -AfT'expr- 
n 


(3.2.3) 


proof:  Choose  /- n  and  y-f„  in  (3.2.1),  and  note  that  Z~(„  <  0.  □ 
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Corollary  2:  Upper  bound  on  Ao"expr.  If  fo  <  i,  for  each  /  —  0.1 . »,  then 

Ao'exp,  <  —  •l,"~lexpT.  (3.2.4) 

k  n 


(3.2.3)  leads  directly  to  a  bound  on  the  relative  error  growth  in  one  step  of  the  standard 
divided  difference  algorithm  (2.1.6).  From  (2.4.5)  the  error  in  ./?(A<j'expr)  relative  to  ld'expT  is 


<o 


Ad‘~lexpT(el,'~1-  cq1-1) 

A<i'expT-(f„-£o) 


The  factor 


ro(r,x) 


W~'ex  pT 
A<i'expT|£„-f0|  ' 


(3.2.5) 


which  we  call  the  growth  factor,  controls  the  growth  of  the  relative  error  in  computing  Ad'expT 
by  (2.1.6).  Clearly  when  ro(r;jc)  is  small,  the  relative  error  growth  is  small;  conversely  when  it 
is  large,  the  relative  error  growth  may  be  large.  By  (3.2.3)  the  growth  factor  satisfies 


rg(T;x)  ^ 


n 

T(f.-fo) 


(3.2.6) 


when  ^  f ,  for  all  /. 

(3.2.6)  illustrates  the  dependence  of  error  growth  on  the  three  factors  mentioned  at  the 
beginning  of  this  section.  It  also  permits  us  to  bound  propagated  errors  in  the  divided 
difference  table  in  the  manner  shown  in  Fig.  2.2.3.  We  illustrate  this  in  Fig.  3.2.1  for  a  single 
initial  relative  error  e  in  \(jexpTS  Each  element  of  the  bottom  diagonal  satisfies 

l«d'l  <  T-^DleHn^-lo)]-1. 

j- i 

Note  that  for  equispaced  abscissae,  say  f, - yfi,  we  have  |«d'|  <  |e|/(rS)",  and  the  factorial  can- 
cells. 


tThe  abscissae  are  •  •  •  <  f_j  <  f_2  <  £-i  <  fo  <  f  i  <  #2  <  #3  <  £«  •  •  • ,  in  contrast 

with  our  usual  numbering. 


.  » 

s 


*-4 

•2, 

«2, 

<2: 

«2, 

*23 

<-i 

,  2 

4 

€ 

<o 

«-i 

€2, 

«-2 

€2j 

<0 

3 

*2. 

_ _ n 

«0 

«n 

1 

W  «4t 

t 

Theorem  3:  Suppose  /3  <  <  y  for  all  abscissae  0  ^  n.  Then  there  exist  values 

f €  1/3,  y],  0  <  j  ^  /?,  such  that  for  any  value  of  f  and  every  index  i 

isr'exp,  -«-«,+  f  *  i  f-V-r  l3'2'7> 


>  i  . 


proof:  By  the  chain  rule  for  differentiation,  (3.1.1),  and  (2.1.10), 


-A"exp(T(x-fw)l  -  2*{(f,-fM"'Hexp(T(jir-£«).  r(f,-f)]} 


-  e-T{r-(',+l,X{(«/-f)-A',+,expt(.v.^)) . 

t-0 


Combining  this  with  (3.2.2)  yields 

l.-r'exp,-  (■i  +  «-f,)-A^xpT  +  -*-£l(f rfl’i'^Uf,)!. 
T  T  /-0 

Now  by  Theorem  2,  for  each  j  there  exists  a  f,  €  (/3,  y]  such  that 


'  '***&&& 
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A"+'exprU*/> 


«+l  +  T(f,-f/) 


Ao'expr.  □ 


Setting  /  -  n  and  f  f «  £„  for  all  0  <  j  ^  n  yields  a  sharper  inequality  than  (3.2.3)  when 
f  „  is  the  largest  abscissa. 


Corollary:  If  ^  for  each  /  -0, 1 . w,  then 

I 


ii- 1 


,To  n+l+T<f 


Ao'exp,.  3*  — — r-Ao"  expr. 
w  + 1 


(3.2.8) 


Even  better  inequalities  can  be  derived,  but  at  the  sacrifice  of  simplicity.  We  also  note  that 

«  j 

(?o'»+l+T<f/-f/>  *  1 

because  the  left-hand  side  of  (3.2.7)  and  Ao'exp,  are  independent  of  f,  thus  giving  a  relation 
amongst  the  f ). 


example:  For  evenly  spaced  data  points  real  exponential  divided  differences  can  be  presented 

analytically.  Let  X  —  {£o>fo+28,£o+48 . f0+2/r8l.  where  25  is  the  spacing. 

Then 


Ao’expr 


1  1  cr«0*.,a>^sinh(r5) 

n'.  26  n\  5 


(3.2.9) 


This  expression  yields  very  accurate  divided  differences,  especially  if  we  have  avail¬ 
able  a  good  routine  to  evaluate  the  function  Sh(f)  =  sinh(f)/£.  Fig.  3.2.2  compares 

divided  differences  Ao'expT  for  rr-0, 1 . 24  computed  according  to  (3.2.9)  and 

the  standard  algorithm  (2.1.6).  r  —  1,  £o-*0  and  28**  1.  The  initial  values  Anexp 
are  rounded  to  seven  digits,  and  all  arithmetic  is  performed  in  greater  precision  in 
order  to  isolate  error  growth  due  to  initial  errors. 

The  table  in  Fig.  3.2.3  compares  the  actual  error  growth  per  step  in  using  the  stan¬ 
dard  scheme  with  bounds  derived  from 


Uol  <  +  rg(r,x-)-(|«,"“'l  +  kr'i) 


under  the  assumption  |«r~' |  -  |«o |.  That  is,  l«ol  < where 


P  =  1  +2/-o'(t;.v). 
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/? 

Ad'exp  by  (3.2.9) 

Ao'exp  by  (2.1.6) 

Relative  error  in 
An  exp  by  (2.1.6) 

0 

1.000000 

1.000000 

0.0 

1 

1.718282 

1.718282 

0.0 

2 

1.476247 

1 .476246 

6.774E-7 

3 

8.455359E-1 

8.455363E-1 

-4.731E-7 

4 

3.632173E-1 

3.632166E-1 

1.927E-6 

5 

1.248219E-1 

1.248225E-1 

-4.807E-6 

6 

3.574655E-2 

3.57461  IE-2 

1.231E-5 

7 

8.774665E-3 

8.77481  IE-3 

-1.664E-5 

S 

1.884669E-3 

1.884642E-3 

1.433E-5 

9 

3.598214E-4 

3.5981 86E-4 

7.782E-6 

10 

6.182746E-5 

6.1831 18E-5 

-6.017E-5 

11 

9.657909E-6 

9.655845E-6 

2.137E-4 

12 

1.382918E-6 

1.383742E-6 

-5.958E-4 

13 

1.827879E-7 

1.825252E-7 

1.437E-3 

14 

2.243437E-8 

2.249849E-8 

-2.858E-3 

15 

2.569905E-9 

2.558588E-9 

4.404E-3 

16 

2.759888E-10 

2.772152E-10 

-4.444E-3 

17 

2.789568E-11 

2.787636E-1 1 

6.926E-4 

18 

2.662925E-12 

2.645633E-12 

6.494E-3 

19 

2.408240E-I3 

2.422432E-1 3 

-5.893E-3 

20 

2.069018E-14 

2.148605E-14 

-3.847E-2 

21 

l  .69293 1  E-l  5 

1.349106E-15  * 

2.031  E-l 

22 

1.322242E-16 

2. 1 1 9468E- 1 6 

-6.029E-1 

23 

9.878198E-18 

-3.198523E-18 

1.324 

24 

7.072304E-19 

2.272222E-18 

-2.213 

Fig.  3.2.2:  Ao'exp  computed  from  initial  values  rounded  to  7  digits. 


*  ■  • 
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Error  growth  factors 

P 

Average  error  growth  per  step 

=  2 

Error  growth  bound  using  rg(r\x) 

!+2/(e'-l)-2.!6 

Error  growth  bound  using  (3.2.8)  to  bound  r/j(r: x) 

<  1  +  2!og2  —  2.4 

Error  growth  bound  using  (3.2.3)  to  bound  r/j(r:x) 

3 

Fig.  3.2.3:  Relative  error  growth  and  bounds  for  divided  differences  in  Fig.  3.2.2. 
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Algorithm:  Taylor  series  algorithm  for  AoexPr- 


1.  Initialize  AoI<i  "  1  -  Aq's*  ”  ~~~ .  for  each  k  -0. 1. 


2.  Fory-0, 1,2,  ... 

Ao‘tr*+l  -  (^-a)-A^tr*  +  A<r'u+‘ 


fiTaTj+k  +  \ 

A0\+*+i  -  AoV*  +  for  each  *-0.1.  . 
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£****&?$  t*"? 


The  algorithm  produces  ail  the  values  Adexp,,  0  <  k  <  n,  as  a  bonus  because  the  quantities 
needed  to  form  them  are  intermediate  in  forming  just  lo'expr.  Also,  the  coefficient  evaluations 
can  be  performed  iteratively  for  greater  efficiency. 

We  now  wish  to  examine  how  the  error  in  computing  lo"expT  by  this  algorithm  is  affected 
by  the  parameter  r,  the  order  /»,  and  the  spread  in  the  data  points  —  fn.  Since  lower 
order  divided  differences  play  no  role  in  the  series  computation,  the  subject  of  propagating  ini¬ 
tial  errors  is  not  relevant,  instead,  we  examine  the  effects  of  round-off  errors  in  each  step  of 
the  series  computation  and  obtain  an  overall  error  bound. 

The  algorithm  involves  many  inner  product  computations.  We  consider  two  possible  error 

n 

conditions.  In  the  first,  the  computed  inner  product  _/?2(£a,/3,)  satisfies 

,-0 


[fli -  Zo-0,1  <  «2>kl 

/— 0  »“0  *“0 


(3.3.3) 


for  all  n.  The  error  analysis  here  is  based  upon  methods  presented  by  Wilkinson  [1 963 1  who 
takes  c  as  1.06  times  the  machine  precision.  The  error  bound  (3.3.3)  holds,  for  example,  when 
all  additions  are  performed  in  double  precision  arithmetic  (hence  the  subscript  2  in  fi:)  and 
rounding  to  single  precision  is  done  only  when  the  summation  is  completed. 

n 

In  the  second  case,  the  computed  sum  /?(£a,/3,)  satisfies 

,-o 

[/K^cttP i)  -  2*0,18,1  <  «£(/i  +  2-f)|a,0,|.  (3.3.4) 

—0  ,-0  1-0 

This  bound  holds  when  the  entire  summation  is  performed  in  single  precision  arithmetic. 

The  following  error  bounds  for  Adexpr  are  established  in  Appendix  B.  Under  the  double 
precision  condition  (3.3.3) 


i/72U<j'expr)  -  Ad'exprl  <  e(2  +  r0/2)er#/7-1-^- 


(3.3.5) 


and  under  the  single  precision  condition  (3.3.4) 


L/WAo'exPr)  -  Ao'exPrl  <  «(m  +  n  +  7  +  T9/2)eT,in^-~ 


(3.3.6) 


The  factor  r"eTa/nl  is  l"expT  evaluated  for  data  points  confluent  at  a.  m+ 1  is  the  number  of 


i 


aggggiww 
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terms  added  in  our  computed  partial  sum  of  the  Taylor  expansion  and  is  chosen  so  that,  say, 
llo'exp,  -  A<S'sw+J  -  |er“  £  -rf— T^tr'l  <  . 


m  depends  on  c,  r,  and  0  in  a  complicated  way;  only  a  general  estimate  can  be  given.  For 
example  Appendix  B  shows  that  when  «-  10-U  and  r0  <  2,  m  >  16  satisfies  the  above  condi¬ 
tion. 


The  bounds  (3.3.5)  and  (3.3.6)  are  converted  to  relative  error  bounds  by 

<  A°exp'  < 


which  follows  from  A"expr  being  increasing  in  its  abscissae.  Then  because  <  <*  <  f 


r"eT' 

n[ 


<  er<a-<°,-A(5'expT  -  eT*/2'Ao'expT. 


Relative  error  bounds  for  the  Taylor  series  algorithm.  The  relative  error  in  represent¬ 
ing  Ao'expT  by  its  computed  value  satisfies 

I tol  <  e( 2  +  T0/2)er"  (3.3.7) 

for  double  precision  accumulation  (3.3.3),  and 

|«o'J  <  «(m  +  n  +  7  +  T0/2)er*  (3.3.8) 

for  single  precision  accumulation  (3.3.4). 


The  relative  error  bound  in  the  first  case  does  not  depend  on  n.  In  both  cases  it  is 
increasing  in  the  "spread"  t9.  These  error  bounds  are  uniform  in  the  sense  that  if  the  Taylor 
series  algorithm  were  used  to  compute  any  other  divided  difference  of  the  table  (any  A,‘expr 

for  A  —  0. 1 . n  and  /  —  0. 1 . n-k ),  a  smaller  error  bound  would  result.  This  follows 

from  the  ordering  condition  (o^(,  ^  Error  bounds  for  A,‘expr  would  either  involve 
replacing  n  by  k  ^  n  or  9  by  a  smaller  number. 

example:  In  Fig.  3.3.1,  8-th  order  divided  differences  correct  to  7  digits  are  given  initially  for 
the  standard  scheme;  the  scheme  is  used  only  to  compute  the  remaining  higher 
order  differences.  The  relative  error  increases  by  a  factor  of  about  3  per  step.  Thus 

l«fl  =  3"“»« 
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n 

abscissae 

f. 

lo'exp. 
correct  to  7 
decimal  digits 

lo'exp  using 
standard  scheme 
after  n  -  8 

lo'exp  using 
Taylor  series 
Algorithm 

lo'exp  using 
Algorithm  of 
§3.4 

0 

-13.0 

2.260329E-06 

2.260332E-06 

1 

-12.5 

2.932648E-06 

2.917451  E-06 

2.932650E-06 

2 

-12.0 

1.902471  E-06 

1.905751  E-06 

1.902472E-06 

3 

-11.5 

8.227822E-07 

8.215427E-07 

8.227824E-07 

4 

-11.0 

2.668782E-07 

2.669856E-07 

2.668782E-07 

5 

-10.5 

6.925 181 E-08 

6.925365E-08 

6.925 180E-08 

6 

-10.0 

1.497504E-08 

1.496362E-08 

1.497505  E-08 

7 

-9.5 

2.775608E-09 

2.777049E-09 

2.775609E-09 

8 

-9.0 

4.501 490E- 10 

4.501490E-10 

4.501040E-10 

4.50 1 49 1 E- 1 0 

9 

-8.5 

6.489361  E-tl 

6.489360E-1 1 

6.490737E-11 

6.489364E-1 1 

10 

-8.0 

8.419572E-12 

8.419580E-12 

8.420343E-12 

8.41 9573E- 1 2 

11 

-7.5 

9.930829E-13 

9.930800E-13 

9.930225E-13 

9.930829E-13 

12 

-7.0 

1.073723E-13 

1.073727E-13 

1.073735E-13 

1.073724E-13 

13 

-6.5 

1.07161  IE-14 

1.071609E-14 

1.071603E-14 

1 .071 6 1 1 E- 1 4 

14 

-6.0 

9.93 1098E- 16 

9.93 1 271 E-16 

9.930869E-16 

9.931 100E-16 

15 

-5.5 

8.590019E-17 

8.590612E-17 

8.590039E-17 

8.590021  E- 17 

16 

-5.0 

6.965660E-18 

6.951898E-18 

6.96561 2E- 1 8 

6.965660E-18 

17 

-4.5 

5.316202E-19 

5.402473E-19 

5.316202E-19 

5.31 6201 E- 19 

18 

-4.0 

3.831926E-20 

3.486964E-20 

3.831920E-20 

3.831926E-20 

19 

-3.5 

2.616686E-21 

3.650387E-21 

2.616686E-21 

2.616687E-21 

20 

-3.0 

1.697500E-22 

-6.986900E-23 

1.697499E-22 

1.697500E-22 

21 

-2.5 

1.048766E-23 

5.010054E-23 

1.048766E-23 

1.048766E-23 

22 

-2.0 

6.185062E-25 

-1.792933E-24 

6. 1 85061 E-25 

6.185063E-25 

23 

-1.5 

3.489027E-26 

- 1 .23641 9E-24 

3.489027E-26 

3.489027E-26 

24 

-1.0 

1 .8861 72E-27 

6.496526E-25 

1 .886 1 72E-27 

1.886172E-27 

25 

-0.5 

9.788799E-29 

-1.931645E-25 

9.788798E-29 

9.788797E-29 

Fig.  3.3.1:  Top  row  of  lexp  using  several  methods. 


where  e^5x  10-7,  whereas  the  growth  factor  bound 

1  a.  h 
P~  + 

discussed  along  with  Fig.  3.2.3  gives  a  bound  of  5  for  the  increase.  The  Taylor 
scheme  yields  good  results  for  a -25  because  the  data  points  are  symmetrically 
placed  about  the  expansion  point  a  -  —7;  however,  the  lower  order  differences  have 
less  relative  accuracy  than  l<J5exp.  The  relative  error  bound  (3.3.7)  with  9-12.5  is 


(ell  <  2.3  x  I06e . 


For  the  lower  order  differences  (small  n),  this  overestimates  |<d'i  by  a  factor  of 
about  10.  Without  a  correct  value  of  the  divided  difference  to  compare  with, 
bounds  such  as  the  above  must  be  accepted  as  the  uncertainty  in  the  computed 
divided  differences,  for  ail  n. 

The  example  shows  that  both  the  Taylor  series  algorithm  and  the  standard  scheme  (even 
with  some  low  order  differences  initially  provided)  may  produce  Ad'expr  with  large  relative 
errors  when  ri i  is  neither  large  nor  small.  The  algorithm  presented  in  the  next  section  is 
designed  to  deal  with  this  intermediate  situation. 
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3.4  A  scaling  and  squaring  method. 

At  the  end  of  the  last  section  we  saw  that  there  are  situations  where  neither  the  standard 
scheme  nor  the  Taylor  series  algorithm  yields  a  value  of  Ao'exp,  with  small  relative  error  for  all 
n  of  interest.  We  present  here  a  third  approach  for  computing  divided  differences  of  the 
exponential  which,  in  many  cases,  gives  significantly  better  error  bounds.  The  method  is  based 
on  the  matrix  function  theorem  for  divided  difference  tables  (§2.6)  and  is  suggested  in  formula 
(2.6.8). 

The  entire  divided  difference  table  is  representable  as  a  function  of  the  special  "step 
matrix"  Z  (2.6.2).  Specifically  for  / -expT, 

Aexp,  -  exp(rZ)  =  eTl  (3.4.1) 

where  the  diagonal  of  Z  consists  of  the  abscissae  f0.£i . €»-  Special  properties  of  the 

exponential  function  are  reflected  in  the  divided  difference  table,  denoted  by  Aexp..  In  particu¬ 
lar  for  any  non-negative  integer  j, 

Aexpr  -  [exp(2~yrZ)P  -  Uexp2_,rP.  (3.4.2) 

a  formula  for  scaling  and  squaring  the  divided  difference  table.  Ward  [1977]  has  suggested  scal¬ 
ing  and  squaring  as  a  method  for  computing  the  exponential  of  a  full  matrix,  whereas  we  pro¬ 
pose  to  use  it  only  in  connection  with  Z. 


example:  With  abscissae  (0. 1 . 2, 3. 4|, 


i 

i  ■. 


1.0000  6.4872E-1 

2.1042E-1 

4.5501  E-2 

7.3794E-3 

1.6487 

1.0696 

3.4692E-1 

7.5019E-2 

■l 

Aexp,/,  - 

2.7183 

1.7634 

5.7198E-1 

4.4817 

2.9074 

t 

1 

7.3891 

to  five  digits.  Squaring  this  matrix 

yields 

1.0000  1.7183 

1.4763 

8.4553E-1 

3.6322E-1  1 

2.2984 

10.908 

34.513 

54.599 


(Aexp./,)2 


2.7182  4.6709  4.0129 

7.3892  12.696 

20.086 


For  example. 


Aoexp  -  £A0*exp.,,-A4  *exp,/. 

A-0 

-  (1.0000)  X  (0.0073794)  +  (0.64872)  x  (0.075019)  +  (0.21042)  x  (0.57198) 

+  (0.045501)  x  (2.9074)  +  (0.0073794)  x  (7.3891) 

-  0.36322; 

and  we  check  by  (3.2.9) 

A04exp  -  -}-(<?'  -I)4  -  -}r  *  (1.7183)4  -  0.36323  . 

4!  4! 

Basic  scaling  and  squaring  algorithm.  For  (3.4.2)  to  be  the  basis  of  a  useful  algorithm,  we 
first  must  obtain  an  initial  divided  difference  table  Aexp2_,T.  In  the  last  section  we  saw  that  the 
relative  error  in  the  Taylor  series  algorithm  decreases  with  the  parameter  r.  By  choosing  j  large 
enough  we  can  make,  say,  the  error  bound  (3.3.7) 

|<tf|  «(2  +  2~('+l>r0)e2~'T'1 

small  for  any  spread  9  in  the  abscissae.  Let  Bj  represent  one  of  the  coefficients  of  e  in  (3.3.7) 
or  (3.3.8),  that  is 

J3;  -  (2  +  2-(/-h'>T0)e2~'r''  or  0>-(«  +  #r  +  7  +  2 (3.4.3) 

The  relative  error  bound  0,«  is  uniform  for  every  element  of  the  divided  difference  table 
Aexp,_,r  computed  by  the  Taylor  algorithm.  Thus 

[/?(Aexp2_,T)  -  AexprJ  <  «0,-Aexpr,T .  (3.4.4) 

The  inequality  holds  element  by  element. r 

tFor  a  matrix  B,  |fl|  denotes  the  matrix  all  whose  elements  are  the  absolute  values  of 
the  elements  of  B,  i.e.  |B|,,/- |B,,|.  Our  notation  B^C  means  that  B,,  <  C,.;  for 
every  i  and  j. 


j  t  » 
I 


In  §2.9  we  remarked  that  it  is  not  necessary  to  compute  an  entire  divided  difference  table 
by  the  series  algorithm.  The  Taylor  algorithm  need  only  produce  the  top  row  of  Aexpr  ,r.  The 

backfill  formula  (2.9.1)  generates  the  remainder  of  Aexprv  Specifically  when  A(?expr,r . 

Ao'cxPj-,,  are  generated  by  the  Taylor  algorithm, 

A ,‘expr,T  -  ({(+*  ,‘_Vexp2.,T  +  A  *_,expr,T  (3.4.5) 

for  fc  —  0.1 . n-i  are  obtainable  successively  for  /  —  1 . 2 . n.  We  show  that  this 

modification  does  not  increase  the  error  bound  (3.4.4). 

example: 


A  exp./. 


6.4872E-1 

2.1042E-1 

4.5501E-2 

7.3794E-3 

1.6487 

1.0696 

3.4692E-1 

7.501 9E-2 

2.7183 

1.7634 

5.7198E-1 

4.4817 

2.9074 

7.3890 

The  top  row  is  from  the  matrix  in  the  previous  example.  The  remainder  of  the  table 
was  filled  in  by  (3.4.5).  For  example. 

A  iJexp./,  —  (4-0),Aoexp./,  +  Aoexpv, 

-  4  x  0.0073794  +  0.045501 

-0.075019 

Lemma:  The  relative  error  bound  on  every  element  of  the  divided  difference  table  is  no 
greater  than  the  largest  error  bound  on  the  elements  in  the  top  row  of  the  table.  This  error 
bound  is  not  increased  when  the  table  is  filled  out  by  the  backfill  scheme  (3.4.5). 

proof:  Consider  one  step  of  the  backfill  routine 

A,"-'exprv  -  (fff-#o>'‘M'expr,T  +  A<T'exp2_,r. 
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Let  «n  be  the  relative  error  in  the  computed  value  of  A<i'expr,r.  The  propagated  error  in 
Ar'expr,T  is 

«r'  A,"-,ex Pj_)r  -  €(i,(|/1-§ o)-A6'expr,r  +  eS~'^n~'ex pr,r. 

Both  |«o I  and  |«o_l|  <6)3/,  by  the  uniformity  of  the  bound  /3,  for  the  Taylor  series 
method.  Thus 

|«r~l|-A|"'~lexpr,T  <  </3,|(€„-fo)-Ad'expr,r  +  A<r'expr,TL 
so 

kr'l  <  *fij. 

When  the  abscissae  are  ordered  <  •  •  •  <  the  sum  in  (3.4.5)  involves  non-negative 

numbers  only  and  the  above  argument  may  be  repeated.  It  shows  that,  considering  only  pro¬ 
pagated  errors,  the  uniform  error  bound  (3.4.4)  holds  also  when  only  the  top  row  of  lexpr,T  is 
computed  by  the  Taylor  series  algorithm  and  the  remainder  of  the  table  backfilled  according  to 
(3.4.5).  □ 

The  outline  of  a  new  approach  for  computing  the  divided  difference  table  Aexpr  is  sum¬ 
marized  as  follows. 


'  i 

' ) 

Error  bounds  and  selection  of  a  scaling  parameter.  An  error  analysis  of  the  algorithm  shows 
how  to  select  j.  The  elements  of  the  divided  difference  table  lexpT  are  non-negative  for  all 
t  >  0.  When  inner  products  involved  in  matrix  squaring  are  accumulated  according  to  the  error 
condition  (3.3.3),  not  depending  upon  n,  we  have  an  element  by  element  error  bound 

\fl2(B2)  -  B2 1  s;  <B2 


Algorithm  1:  Scaling  and  squaring  algorithm  for  Aexpr. 

With  the  abscissae  ordered  such  that  |0  4  §i  <  •  •  •  <  f 

1.  Choose  j  and  form  Aoexp2-,r  f°r  k-  0,1 . n  by  the  Taylor  series  algorithm 

(§3.3); 

2.  Backfill  the  remainder  of  the  table  Aexp2_/r  according  to  (3.4.5); 

3.  Square  the  divided  difference  table  matrix  j  times. 


(3.4.6) 


for  any  non-negative  matrix  S.  such  as  8  -  Jkexp,.v  When  a  computed  matrix  ,/M£>  satisfies 

MB)  -  B  =  E  where  \E\  <  cfi  B .  <3.4.7t 

as  in  (3.4.4),  then 

!(/M«)]-  -  B: )  -  | BE  +  EB  +  £2!  ^  (2c/3,  +  t20;)  B2 . 

Thus  squaring  a  computed  matrix  MB)  yields 

\MMB)\2  -  B2\  «:  \MMB)]2  -  [MB))2\  +  ll/MS)]2  -  b2 ! 
e[/t2(B)}2  +  (2tp,+e:j32)B2 
<  e [(1  +  2/3,)  +  «(2  +  /3,)/3,  + 

«  is  so  small  that  terms  in  e;  are  negligible  when  compared  with  terms  linear  in  e.  We  take 

I/! hWB)]2-  B2\  ^  e(l  +  20,)  B2 .  (3.4.8) 

In  (3.4.8)  S«Aexp,,,r,  so 

«  /72l/f;(Aexp,_,r)l;  =  /MAexp,.  . 

The  first  computed  matrix  square  satisfies 

;./MAexp;.  ....)  -  Aexpr,,.,,r|  ^  e(l  +  2/3,)-Aexp,-,,.,,r 

-  «l2(j8,  +  1)  -  l]  Aexp,-,,-(,r . 

This  inequality  is  the  same  as  (3.4.7).  but  now  with  B  —  Aexp, and  0  replaced  by 
/J,,_i  =  2(/J  t!>  -1.  Hence  iteratively 

l/MAexpr, -  Aexp2.,,_;,r|  ^  «(4(/3,  + 1)  -  Ij-Aexp,. 

and  after  j  steps 

i.//j(lexpT)  -  AexpT|  <  e[2'(/3,  + D- U-Aexp. .  (3.4.4) 

This  last  inequality  is  a  relative  error  bound  on  the  divided  difference  table  computed  according 
to  the  scaling  and  squaring  algorithm.  The  bound  depends  on  the  Taylor  series  bound  through 
0  and  increases  exponentially  in  j,  the  number  of  squaring  operations  performed. 
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It  is  clear  both  from  the  error  bound  (3.4.9)  and  the  amount  of  work  involved  in  squaring 
matrices  that  we  wish  to  choose  j  as  small  as  possible.  However  0,  increases  as  /  decreases. 
We  demonstrate  here  how  to  select  j  to  minimize  the  bound  in  (3.4.9).  With  j  so  chosen,  we 
obtain  an  expression  of  the  form  ktQ  for  the  relative  error  bound.  The  value  of  the  constant  k 
depends  on  the  specific  assumptions  made  in  bounding  round-off  errors  in  the  Taylor  algorithm 
and  the  matrix  squaring. 

We  want  to  minimize  the  coefficient  in  (3.4.9);  that  is  we  want  to  choose  /  to  minimize 
the  expression 


g(j)  =  2;(/3,  + 1)  -  1  -  2,l(2  +  2~,'+uTtf)er'r"-t- 1]  -  1  . 

Define  cr  s  2"'t0.  r  and  9  are  fixed  here,  so  minimizing  g(j)  in  /  is  equivalent  to  minimizing 

g(tr)  S  —  ((2  +  o-/2)e"+ ll 
cr 

in  o’.  The  minimum  is  g(<r)  =  7.7885  which  occurs  for  cr  =  0.9606.  2~'tO  —  0.9606  is  probably 
not  true  for  integer  values  of  j.  Nevertheless  for  integer  j,  the  o-  -  2~'t9  yielding  the  smallest 
value  of  g(tr)  must  satisfy  cr0  <  2"'r0  <  2tr0  where  g(cr„)  —  ^  (2cr0) . 


Fig.  3.4.1:  Graph  of  g(<r)  showing  (cr0, 2cr0]  is  the  largest  interval  containing  ct=~T't9 

for  just  one  integer  value  of  j  . 

o-o  =  0.6646,  so  0.6646  <  2~'t9  <  1.3292.  The  minimizing  j  is  the  smallest  non-negative 
integer  satisfying 

2~’r9  <  1.3292.  (3.4.10) 

For  all  <r  €  (0.6646,  1.3292).  g(cr )  <  4U  -3292).  We  are  assured  that  for  the  above  choice  of  j. 
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£(2“'rtf)  <  £(1.3292)  ~8.3259  =  k2.  Then 

g(j)  -  t9 -gicr)  —  1  <  k 2rtf  ~  1  . 

Hence 

|,//2(AexpT)  -  AexpT|  ^  €[K2rtf-  1]  lexpr  (3.4.11) 

when  the  scaling  parameter  j  is  chosen  according  to  (3.4.10).  This  bound  may  fail  when  rtf  is 
very  small;  in  this  case  j  -0  and  the  Taylor  series  bound  (3.3.7)  is  appropriate. 

By  a  similar  argument  we  derive  a  bound  like  (3.4.1 1)  for  the  single  precision  error  condi¬ 
tion  (3.3.4).  In  this  case  we  write 

(3,  -  (fl  +  n-y,  S  (,„  +  *  + 7 +  2-,'+"rtf>c'r'r,\ 

which  is  consistent  with  the  other  Taylor  series  error  bound  (3.3.8)  The  error  in  matrix  squar¬ 
ing  satisfies 

IjHB1)  -  B2 1  ^  *(n+\)B2 
for  any  non-negative  matrix  B.  From  (3.4.7), 

MB)  -  B  =  E  where  )£)  e(/i  +  l)-y,fl  . 

Replacing  e  by  «(n  +  l)  and  /3 ,  by  y ,  in  our  arguments  leading  up  to  (3.4.9), 

L/?(lexpr)  -  JtexpT|  <  «(//+l)l2'(y,  +  1)  -  1]-Aexpr.  (3.4.12) 


In  Appendix  B  we  show  that  when  «  ^  10-7,  m  can  be  taken  as  small  as  10.  We  assume 
also  all  first  order  divided  differences  are  computed  by  a  special  formula  as  in  (3.1.2).  so  our 
bounds  here  are  applied  only  when  n  ^  2.  Also,  j  will  be  such  that  2 ~"*ut9  <  1.  Hence 


?/*■(!  +  • 


+  6  +  2_0+"rtf 


)e2"™  <  le1'" 


(3.4.131 


As  before  we  want  an  integer  j  to  minimize  the  expression 


2'(7e2''T*+  1)  -  l  . 


This  is  minimized  when  j  is  the  smallest  non-negative  integer  satisfying 


•  ~  .  '  v 
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53.4 

2 ~'t#  <  1.4542. 

For  this  j 

L/7(Aexpr)  -  Aexpr|  <  ) [k|T«  —  llAexp, 

where  k t  “21.2950.  The  following  box  summarizes  these  bounds. 


(3.4.14) 


(3.4.15) 


Scaling  and  squaring  error  bounds.  For  double  precision  accumulation  of  inner  products 

(3.3.3), 

|y?2(dkexpr)  -  Aexpr|  <  e[x2r0- l]Aexpr 

where  x2  —  8.3259,  9  is  the  maximum  spread  in  the  abscissae,  and  the  number  of  squarings 
j  is  the  smallest  non-negative  integer  such  that 

2_,r 9  1.3292. 

For  single  precision  accumulation  (3.3.4), 

L/7(Aexpr)  -  Aexpr|  <  e(n+l)[ieiT0-  1]  Aexpr 
where  «|  -  21.2950  and  j  is  the  smallest  non-negative  integer  such  that 

2~>t 9  <  1.4542. 


example:  The  entries  in  the  right  hand  column  of  Fig.  3.3.1  were  computed  by  scaling  and 
squaring.  The  bound  (3.4.11)  has  the  coefficient 

k2t9-  1  -  8.3259  x  1  x  12.5-1  =  103, 

and  logiol03  =  2.01.  This  indicates  a  loss  of  two  decimal  digits,  at  most,  in  all  the 
divided  differences  computed. 


Modify...  scaling  and  squaring  algorithm.  The  algorithm  can  be  made  more  efficient  by 
extending  our  use  of  the  backfill  scheme.  Squaring  a  (n+1)  x  («+ 1)  triangular  matrix  (such  as 
Aexp,.,,.)  requires  (/?+3)(/j+2)(rH-l)/6  multiplications.  The  j  squarings  needed  to  get  Aexpr 
from  Aexpr,r  involve  0(jnll(>)  operations.  This  operation  count  can  be  reduced  to  0(  jn:). 

Once  the  top  (0-th)  row  of 772(Aexp2_(,_,,r)  is  computed  from  squaring  y?2(Aexp,_,r),  the 
backfill  scheme  (2.9.1)  will  generate  the  remainder  of  y/2(<Aexp2_,,_i,T)  in  exactly  the  same 
manner  we  generated  the  remainder  of./?2(Aexpri).)  given  its  top  row  by  the  Taylor  series  algo- 
•ithm  Because  relative  errors  in  the  elements  in  this  top  row  of  y?2<Aexp2_a,_i,r)  are  uniformly 
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bounded  by  e(l  +  2/3,), 

L//2(Aexpr„-,,T)  -  Aexp2-„-nT|  ^  «(1 +2/8,)-Aexp2-„-i,r 

for  the  entire  table.  This  idea  holds  for  all  the  squarings.  We  notice  also  that  the  uniform 
bounds  are  exactly  those  employed  in  the  analysis  of  the  scaling  and  squaring  algorithm.  Thus 
all  our  just  derived  error  bounds  are  applicable  when  the  matrix  squaring  is  modified  in  the 
above  manner.  The  same  argument  holds  for  single  precision  accumulation  when  e  is  replaced 
by  e(/;  +  l). 

Now,  how  does  the  operation  count  change?  Obtaining  the  top  row  of  a  matrix  square 
requires  (h+2)(w+1)/2  multiplications.  Backfilling  the  rest  of  the  matrix  requires  one  multipli¬ 
cation  per  element,  or  «(r/+l)/2  multiplications.  Thus  computing  each  matrix  square  by  this 
modified  method  requires  (n+1)2  multiplications,  compared  with  (n+3)(«+2)(»+l  )/6  for  the 
direct  squaring  approach.  This  is  an  improvement  for  all  n  >  1. 

Algorithm  2:  Modified  scaling  and  squaring  algorithm  for  Aexpr. 

With  the  abscissae  ordered  such  that  fo  <  |j  ^  •  •  •  < 

1.  Choose  j  according  to  (3.4.10)  or  (3  4.14)  and  form  A0‘exp,-.r  for  it -0.1 . n 

by  the  Taylor  series  algorithm  (§3.3); 

2.  Backfill  the  remainder  of  the  table  Aexpr  ,r;  j 

3.  Square  the  divided  difference  table  matrix  Aexp,^  ,r  j  times  in  the  modified  manner;  J 

compute  the  0-th  row  of  the  matrix  square  and  then  backfill  the  remainder  of  the  | 
table.  I 


3.5  A  hybrid  algorithm  for  the  divided  difference  table  Aexpr. 


We  have  now  presented  three  quite  different  methods  for  computing  A»expr 

=  (A"expr)(£o . £„):  (I)  standard  divided  differences,  (2)  Taylor  series,  and  (3)  scaling 

and  squaring.  These  algorithms  have  complementary  error  propagation  properties,  but  they 
vary  in  computational  efficiency.  We  summarize  here  these  two  aspects  of  each  algorithm  and 
present  a  hybrid  algorithm  which  may  be  used  when  none  of  the  above  alone  is  satisfactory  for 
computing  an  entire  table.  For  our  hybrid  algorithm  we  give  error  bounds  which  depend  only 
on  the  order  of  the  divided  differences  computed:  these  bounds  are  independent  of  the  choice 
of  abscissae  and  parameter  r. 

(I)  Standard.  The  propagated  relative  error  ed'  in  a  typical  step  of  the  standard  algorithm 
satisfies 


*(f|  <  l«r‘l  +  rg(r;x-H|«,"-|l  +  l«<Tl|l. 


where  by  (3.2.6) 


rjj(r;x)  < 


n 

r(f„-fo> 


when  ^  f ,  for  all  /. 

When  the  abscissae  are  ordered,  ^  ^  f „,  and  all  initial  relative  errors  e,n  in 

the  function  values  A,°expr  =  e,f  are  uniformly  bounded  (that  is  |e,°|  for  all  / ),  we  obtain 
a  simple  bound  on  eo'.  Let  <6  represent  the  minimum  separation  of  the  data  points,  that  is 

<t>  =  min  (£,  The  relative  error  in  all  first  order  divided  differences  satisfies 

i  <  - « " 

|«,*|  «<!  +  2/t*>.  /-0.1 . n- 1 


Continuing,  we  obtain 


!«ot  <  «Iln  +2*/^>  <3  5.1 ) 

4-1 

This  bound  is  decreasing  in  r*.  Given  the  (initial)  n+1  exponentials  e  ( .  computing  An'expr. 
in  fact  the  entire  divided  difference  table,  requires  only  «(n+l)/2  divisions. 

(2)  Series.  The  Taylor  series  method  (§3.3)  needs  2/?+2  multiplications  (exclusive  of 
coefficient  evaluations)  to  add  a  new  term  to  the  partial  sums  we  form  for  each  An'exp.. 
k  »0. 1 . n.  When  each  partial  sum  has  m+1  terms,  the  total  is  2m(»»+l)  multiplications 
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(the  initial  term  requires  no  inner  products)  to  get  Ao'exp,,  as  well  as  all  Ao‘expr  for 
A- -0,1 . n. 

(3)  Scaling  and  squaring.  Squaring  a  (n  +  1)  x  (n+1)  divided  difference  table  by  the  special 
method  described  in  §3.4  requires  (n-t-1)2  multiplications,  and  we  do  this  j  times.  For  any  but 
small  »,  the  table  squaring  dominates  the  rest  of  the  calculation. 

Fig.  3.5.1  summarizes  this  information  on  bounds  and  operation  counts.  The  error 
bounds  listed  in  the  second  column  assume  that  inner  products  satisfy  the  double  precision 
error  condition  (3.3.3),  that  is 


!A’(Xa,/3,)  -  ^  «II«Al  • 

i  “0  /  i  "O 

Those  in  the  third  column  reflect  the  single  precision  error  condition  (3.3.4),  namely 
L/7(]Ta,0,)  ~  2>,0,|  <  e£(n  +  2-/)|a,/3,| . 

i  »0  i  “0  i  *-0 

The  bounds  depend  on  the  minimum  separation  <f>  and  the  maximum  spread  9  in  the  abscissae. 
The  constants  k2  and  *i  in  the  bottom  entries  depend  on  the  details  of  the  arithmetic  in  the 
scaling  and  squaring  algorithm.  In  §3.4  we  derived  the  values  x2  —  8.3259,  and  *|-  21.2950 
when  «  — 10-7. 


Method 

Relative  error  bound  coefficients 

Operations 

Double  precision 

Single  precision 

Standard  algorithm  with 

minimum  separation  6 

JJU  +2A/r<£) 

A-l 

lid  +2fc/r0) 

^-i 

-1-* 

Taylor  series,  using  tn+1 
terms,  with  spread  9 

(2  +  T0/2)er* 

(m  +  n  +  l  +  T9/2)er'1 

—  2mn 

Scaling  and  squaring 
with  spread  9 

k2t9  —  1 

(tf  +  l)[K|T0  —  11 

-  jn2 

Fig.  3.5.1:  Summary  of  three  divided  difference  algorithms  for  Ao'exp. . 


Decision  criteria  and  the  hybrid  algorithm.  Our  error  bounds  suggest  a  hybrid  algorithm:  com¬ 
mit*  a  II  divided  differences  having  closely  Clustered  abscissae  bv  seal  ins  and  snuarine.  and  the 
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remainder  of  the  table  by  the  standard  formula.  The  Taylor  series  is  a  special  case  of  scaling 
and  squaring  with  /  —  0.  The  operation  counts  suggest  employing  scaling  and  squaring  on 
divided  differences  of  the  smallest  practical  order.  Our  desire  for  good  accuracy  and  our  desire 
for  efficiency,  however,  are  in  conflict.  Here  we  lean  towards  the  former  in  presenting  criteria 
for  deciding  which  method  to  employ  when  computing  a  particular  difference  in  the  table;  we 
develop  an  overall  error  bound  which  holds  for  any  sequence  of  abscissae  and  parameter  r. 

A  simple  criterion  is  to  use  the  "best"  method  to  compute  each  divided  difference  in  the 
table.  By  best  we  mean  that  method  (either  scaling  and  squaring  or  standard)  which  yields  the 
smallest  relative  error  bound  for  that  particular  divided  difference  we  are  considering.  All  lower 
order  divided  differences  are  assumed  to  have  been  computed  already  by  the  best  method,  or 
by  a  special  formula. 

For  example  in  computing  Ao'expr  with  double  precision  accumulation,  we  use  scaling  and 
squaring  when 


(*<2T0-1)«  <  (l+2/r0)e, 

and  the  standard  algorithm,  otherwise.  Here  is  both  the  spread  and  the  minimum 

separation  in  the  data  points.  The  worst  possible  error  bound  for  this  hybrid,  then,  occurs 
when 

k2t  9  —  1  —  1  +  2/r0 . 


This  has  the  solution 

t9  —  (I  +y/l  +  2k2)/k 2  —  T#|  • 

For  k2- 8.3259  (as  derived  in  §3.4),  rtf  |  ==  0.62  and 

k2t9}  —  1  =  1+  2/rtf  |  =  4.20 . 

Thus  the  relative  error  «d  is  bounded,  j«ol  ^  4.20«,  when  the  "best"  method  is  used.  This 
bound  does  not  depend  on  the  abscissae  or  r,  only  on  the  value  of  k2.  We  have  obtained  a 
bound  independent  of  the  abscissae  and  r  when  we  use  scaling  and  squaring  for  t9  <  t9 ,  and 
the  standard  formula  for  r9  ^  r9{.  This  is  a  simple  criterion  for  deciding  which  method  to 
employ. 
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re*,  e  r»r>  o  r 


Fig.  3.5.2:  Uncertainty  in  computed  values  of  Aoexpr. 


It  is  important  to  note  exactly  what  our  criterion  means.  The  case  tO-t#,  does  not  mean 
that  the  two  methods  are  equally  accurate,  only  that  our  convenient  error  bounds  for  each 
method  are  equal.  Each  bound  may  be  viewed  as  our  maximum  uncertainty  in  the  computed 
A,jexpr  when  the  appropriate  method  is  used.  Thus  when  r0  —  our  uncertainty  is  equalized 
for  the  two  methods,  and  is  maximized  over  all  tO  for  the  hybrid  method.  The  number  4.20«, 
for  example,  represents  our  maximum  uncertainty  in  the  computed  l(jexpT  when  the  "best" 
method  is  used.  More  refined  error  bounds  using  information  about  r  and  the  abscissae  will 
reduce  the  uncertainty,  but  at  the  loss  of  the  simplicity  we  have  here. 

For  A02expr  and  — f0,  the  relative  error  for  the  standard  formula  is  bounded  by 
(k2t9i  —  1 )  ( 1  +4/r0)«.  We  use  scaling  and  squaring  when 

(k2t0  — 1)«  <  (k2t0]  —  1 ) ( 1  +4/r0)e  . 

The  largest  error  bound  occurs  when  equality  holds.  Let  this  happen  for  t0  —  t02,  thus 

KiT0j—  I  —  (*jT0|  —  1 ) ( 1  +  4/t02)  . 

This  procedure  may  be  followed  for  all  divided  differences.  For  Jto'exp,  we  obtain  the 
recurrence  in  t9„, 

kjt9„  —  1  »  (Kjr0„_i  —  I ) ( I  +  2 nlr0„) .  (3.5.2) 

The  criterion  for  scaling  and  squaring  in  computing  4<fexpr  is 

r((„-(o)  <  t0„,  (3.5.3) 


m 


Error  bound 
coefficient 


Bound  on  decimal 
digits  lost 


«(//-  3) 


logio[K2n(/>-3)  —  ll 


89.67 

109.32 

131.00 

154.69 

180.39 

208.11 

237.85 

269.59 

303.35 

339.11 

376.89 

416.68 

458.47 

502.27 

548.08 

807.23 

1,116.50 

1,475.87 

1,885.31 

2,344.82 

3,413.96 

4,683.24 

6,152.62 

7,822.07 

9,691.59 


4.20 

13.71 

33.54 

66.59 

114.57 

178.38 

258.51 

355.19 

468.55 
598.66 

745.55 
909.22 

1,089.68 

1,286.92 

1.500.94 

1,731.73 

1.979.28 
2,243.58 
2,524.63 
2,822.43 

3.136.95 
3,468.21 
3,816.18 

4.180.88 

4.562.29 

6,719.91 

9.294.88 
12,286.96 
15,695.94 
19,521.72 

28,423.30 

38,991.19 

51,225.09 

65,124.80 

80,690.18 


Fig.  3.S.3:  Error  bounds  and  decision  criteria  for  hybrid  algorithm. 

and  the  relative  error  is  bounded  by  (K2r0„-l)e.  Since  *2r0o-l  3 1,  the  recurrence  (3.5.2) 
can  be  evaluated  to  yield  t9„  for  any  n.  These  values,  along  with  corresponding  relative  error 
bounds,  are  listed  in  Fig.  3.5.3  when  *2-  8.3259. 

We  can  now  summarize  a  complete  algorithm  for  computing  divided  differences  of  the 


exponential  with  real  abscissae. 


Sr  rSLS 
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Algorithm:  Hybrid  algorithm  for  Aexpr. 

1.  Compute  for  each  <—0.1 . /». 

2.  For  A  -  1.2 . n  and  /  — 0, 1 . n—k ,  when 

<  r», 

compute  A  ,*expr  by  scaling  and  squaring;  otherwise,  when 

r<§, +*-£,)  >  r»k 

compute  A,Aexpr  by  the  standard  formula. 


The  hybrid  algorithm  requires  us  to  decide  which  divided  difference  scheme  to  use  for 
each  divided  difference.  For  example  in  computing  Ao2expr  by  employing  the  values  in  Fig. 
3.5.3,  when  all  lower  order  divided  differences  have  been  computed  according  to  the  algorithm, 
scaling  and  squaring  is  used  when 

r<«i2-fo)  <  rtfI2=  109.32. 

The  standard  scheme  is  used  otherwise.  The  relative  error  in  our  computed  Ao2expr,  that  is 
«o2,  satisfies 

|€<j2|  <  Oc2t012-1)€  =  909.226. 

log)0(909.22)  =  3  bounds  the  number  of  decimal  digits  lost  in  computing  Ao2expT  by  the  hybrid 
algorithm.  That  is,  when  all  A°expT  are  given  to  10  correct  decimal  digits,  cty,  our  computed 
A<j2expr  contains,  at  least,  7  correct  decimal  digits. 

To  gain  a  better  idea  of  how  the  decision  criterion  t9„,  and  its  associated  error  bound 
*c2r0„-l,  depend  upon  the  order  of  the  divided  difference  <t,  we  bound  solutions  of  the 
recurrence  (3.5.2).  Appendix  C  shows  that 

t9„  <  n1  +  n  +  —  (3.5.4) 

*2 

for  n  ^  l  and  any  *2>0.  Hence  the  relative  error  in  Ao'expT,  computed  according  to  the 
hybrid  algorithm  with  double  precision  accumulation,  satisfies 

|«<{'|  <  (*c2(ffJ  +  n)  +  1)«.  (3.5.5) 

The  relative  error,  then,  increases  in  n ,  at  worst,  as  0(n2).  This  bound  holds  regardless  of  our 
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choice  of  abscissae  and  parameter  r.  For  our  example  with  k2— 8.3259, 

n(n—  4)  <  t0„  <  n(.n- 3) 

when  n  >  17  (Appendix  C).  For  comparison,  the  last  two  columns  of  Fig.  3.5.3  contain  values 
of  /i  (« — 3)  and  logmUj/rf/i— 3)  —  ll.  The  latter  numbers  closely  bound  the  digits  lost  values 
for  large  n. 

Single  precision  decision  criteria.  A  similar  analysis  in  the  case  of  single  precision  accumula¬ 
tion  shows  that  the  hybrid  algorithm  error  bound  behaves  as  0(n}).  Here  the  scaling  and 
squaring  error  satisfies 

|«o'|  <  (n+  1)(*|T0  —  l)e  , 

but  the  bound  on  the  standard  scheme  is  unchanged.  The  same  argument  as  before  leads,  now, 
to  the  recurrence 

(w+1)(k|T0„—  1)  «■  rt(»c |T0„_!  —  1)(1  +2 «/t0„)  .  (3.5.6) 

In  deriving  the  scaling  and  squaring  bound  we  assumed  ail  first  order  differences  are  computed 
by  a  special  method,  therefore  we  require  initially 

2(k|T0)  —  l)  —  I ; 


hence 


is  the  initial  value.  The  table  in  Fig.3.5.4  lists  values  of  the  decision  criterion  r0„  and  its  asso¬ 
ciated  error  bound  for  21.2950,  which  was  derived  in  §3.4. 

We  can  also  show  how  r0„  and  its  associated  error  bound  for  the  hybrid  algorithm  depend 
upon  n  in  the  single  precision  accumulation  case.  From  Appendix  C, 

rd"  *  i"2+(27r_i)w  a5-7> 


for  n  >  1  and  all  *\  >  0;  hence 

|«<$1  <  «(/r-»-l)l-f-K(rr2  +  -  ll 


(3.5.8) 


n 

r0„ 

Error  bound 
coefficient 

Bound  on  decimal 
digits  lost 

n(2n- 5>/3 

j 

logiol(«+n{*r 
•»(2it  — 5)/3  -  1 ))  1 

i 

0.07 

1.00 

0.00 

-1.00 

2 

0.28 

15.11 

1.18 

-0.67 

3 

1.15 

93.95 

1.97 

1.00 

1.91 

4 

3.16 

331.66 

2.52 

4.00 

2.62 

5 

6.58 

835.34 

2.92 

8.33 

3.02 

6 

11.50 

1,707.08 

3.23 

14.00 

3.3 

2 

7 

17.90 

3,041.99 

3.48 

21.00 

3.5 

5 

8 

25.77 

4,930.48 

3.69 

29.00 

3.7 

5 

9 

35.08 

7,460.36 

3.87 

39.00 

3.9 

2 

10 

45.80 

10,717.98 

4.03 

50.00 

7 

1 

11 

57.92 

14,789.00 

4.17 

62.33 

4.2 

12 

71.42 

19,758.68 

4.30 

76.00 

4.3 

2 

13 

86.29 

25,712.06 

4.41 

91.00 

4.4 

3 

14 

102.53 

32,734.11 

4.52 

107.33 

4.5 

3 

15 

120.12 

40,909.77 

4.61 

125.00 

4.6 

3 

16 

139.06 

50,323.94 

4.70 

144.00 

4.72 

17 

159.35 

61,061.56 

4.79 

164.33 

4.8 

18 

180.98 

73,207.55 

4.86 

186.00 

4.88 

19 

203.96 

86,846.88 

4.94 

209.00 

4.95 

20 

228.28 

102,064.51 

5.01 

233.33 

5.02 

21 

253.94 

118,945.43 

5.08 

259.00 

5.08 

22 

280.93 

137,574.66 

5.14 

286.00 

5.15 

23 

309.27 

158,037.22 

5.20 

314.33 

5.21 

24 

338.94 

180,418.14 

5.26 

344.00 

5.26 

25 

369.95 

204,802.47 

5.31 

375.00 

5.32 

1 

30 

545.01 

359,752.88 

5.56 

550.00 

5.5 

6 

35 

753.42 

577,551.81 

5.76 

758.33 

5.7 

6 

40 

995.17 

868,839.39 

5.94 

1,000.00 

5.9 

4  i 

45 

1,270.26 

1,244.258.17 

6.09 

1.275.00 

6.1 

50 

1,578.67 

1,714,452.19 

6.23 

1,583.33 

6.2 

4 

60 

2,295.47 

2,981,746.70 

6.47 

2,300.00 

6.4 

8 

70 

3,145.59 

4,755,889.66 

6.68 

3.150.00 

6.6 

8 

80 

4,129.02 

7,122,053.25 

6.85 

4,133.33 

6.8 

5 

5,245.78 

10,165,412.46 

7.01 

5,250.00 

1 

6,495.85 

13,971,143.98 

7.15 

6,500.00 

7.1 

5 

Fig.  3.5.4:  Single  precision  decision  criteria  and  error  bounds 
for  the  hybrid  algorithm. 


for  the  hybrid  algorithm.  Thus  the  relative  error  is,  at  worst,  0(/i})  in  n.  We  stress  that  this 
bound  holds  for  any  choice  of  abscissae  and  parameter  r.  Further,  in  the  example  with 
*i  -21.2950,  Appendix  C  shows  that 
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2  2  •%  a  ^  2  5 

y/r  -2  n  <  t9„  <  y/r  —  y/i 

for  //  ^  15.  The  rightmost  two  columns  of  Fig.  3.5.4  list  values  of  //(2«— 5)/3  and 
logio((/t+l)U|//(2fl-5)/3  —  1}]. 

The  hybrid  algorithm  demonstrates  that  it  is  possible  to  compute  exponential  divided 
differences  to  a  desired  accuracy.  Our  error  bounds,  particularly  the  digits  lost  bounds,  tell  how 
many  decimal  digits  we  must  carry  in  order  to  be  assured  that  Ao'exp,  has  desired  accuracy.  A 
short  discussion  of  some  useful  modifications  of  the  basic  hybrid  algorithm  follows  in  the  next 
section,  along  with  a  numerical  example  in  which  a  rather  large  divided  difference  table  is  com¬ 
puted. 
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3.6  Comments  on  the  hybrid  divided  difference  algorithm. 

The  hybrid  algorithm  of  the  last  section  demonstrates  that  we  can  compute  high  order 
exponential  divided  differences  with  only  a  modest  loss  of  precision.  For  this  reason,  the  algo¬ 
rithm  is  a  valuable  theoretical  tool.  In  this  section  we  propose  some  modifications  which  make 
its  implementation  more  efficient. 

We  gave  scant  consideration  to  computational  efficiency  when  deriving  the  hybrid  algo¬ 
rithm.  Our  error  bounds  and  decision  criteria  apply  without  reference  to  any  particular 
sequence  of  data  points  or  parameter  r.  As  a  result,  the  algorithm  recomputes  low  order 
divided  differences  when  scaling  and  squaring  is  used  for  differences  whose  "patterns  of  depen¬ 
dence"  overlap.  Also,  the  decision  criteria  are  based  upon  worst  case  arrangements  of  the 
abscissae.  These  arrangements  cannot  be  achieved  since  it  is  impossible  to  arrange  even  three 
points  on  a  line  such  that  their  separations  are  quadratic.  A  relaxed  decision  criterion  may 
greatly  increase  efficiency,  without  sacrificing  accuracy.  We  now  propose  a  possible  modification 
to  the  algorithm  by  introducing  an  arbitrary  criterion  to  cluster  the  data  points. 

Clustering.  Let  g  be  a  positive  increasing  function  of  the  order  k  of  the  divided  difference 
under  consideration.  We  decide  to  use  scaling  and  squaring  to  compute  A,^expr  when 

T<*H.4-f(>  <  *<*>;  (3.6.1) 

otherwise,  we  use  the  standard  formula.  In  addition,  however,  we  do  not  permit  the  computa¬ 
tion  of  overlapping  table  blocks  by  scaling  and  squaring.  For  example,  suppose  the  decision  cri¬ 
terion  (3.6.1)  demands  that  both  A,‘expT  and  A, 'ex p.,  with  /  ^  i+k  <./+/,  be  computed  by 

scaling  and  squaring.  We  compute  only  A/*'_'expr  by  scaling  and  squaring,  regardless  of 
whether  or  not 


r(f,+,-f,)  <  g(j+l-i)  . 

The  picture  in  Fig.  3.6.1  shows  how  overlapping  blocks  may  be  combined.  We  now  speak  of 

the  abscissae  f,.£,+i . f(+,  as  being  "clustered,"  and  refer  to  the  block  of  the  divided 

difference  table  formed  by  A  ,/+/-'expT’s  pattern  of  dependence  as  corresponding  to  this  cluster 
of  abscissae.  So  when  two  clusters  overlap,  they  are  combined. 

When  the  clustering  procedure  is  completed,  the  resulting  clusters  have  no  abscissae  in 
common;  the  corresponding  blocks  in  the  table  do  not  overlap.  This  clustering  depends  on  the 
abscissae,  not  on  the  divided  differences,  and  it  may  be  performed  prior  to  any  divided 
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x,A  •  x,'*1'1 

X 

X  X  x  xj 

X  XXX 
XXX 
X  X 
X 

Fig.  3.6.1:  Clustering  of  overlapping  blocks  in  a  table. 

difference  computations.  Table  blocks  corresponding  to  the  resulting  clusters  are  computed  by 
scaling  and  squaring  with  backfill.  The  picture  in  Fig.  3.6.2  shows  what  the  table  might  then 
resemble.  The  remainder  of  the  table  is  filled  in  by  the  standard  scheme. 

x  x  x  x 
xxx 

X  X 
X 


XXX 
X  X 
X 

XXX 
X  X 
X 

Fig.  3.6.2:  Block  structure  of  a  divided  difference  table  after 
the  scaling  and  squaring  step. 

Our  error  bounds  make  possible  a  quick  a  priori  error  bound  computation.  For  example, 
the  bound  (3.4.11)  may  be  used  when  scaling  and  squaring  is  indicated,  and  the  iterative  bound 
(3.2.10)  when  the  standard  scheme  is  called  for.  We  may  even  wish  to  compute  the  differences 
using  the  decision  criterion  (3.6.1)  and  then  examine  a  postiori  error  bounds.  In  any  event, 
when  these  bounds  are  unacceptable  the  original  hybrid  algorithm  does  guarantee  a  bound  on 
the  error  and  may  be  used  when  more  efficient  alternatives  fail. 


xxx 
x  x 
x 
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Divided  Differences  by  the  Hybrid  Algorithm 

n 

Abscissae 

Clustering 

Ao'exp  by 

Correct  values 

A  priori 

Relative 

algorithm 

to  7  digits 

bound 

error 

0 

-34.5 

0 

1.039538E-15 

1.039538E-15 

0.00 

0.00 

1 

-33.1 

0 

2.268571  H- 1 5 

2.26857 1 E- 1 5 

0.50 

0.00 

2 

32.9 

0 

1.498804E  15 

1.498804E- 15 

2.00 

0.00 

3 

-14.4 

3 

8.015853E-1 1 

8.015853E-1 1 

0.77 

0.00 

4 

-14.4 

3 

6.7551 17E-1 1 

6.7551  1 7E- 11 

0.69 

0.00 

5 

-14.4 

3 

2.879424E-1 1 

2.879424E-1 1 

1.61 

0.21 

6 

-14.4 

3 

8.262803E-12 

8.262803E-12 

1.83 

0.00 

7 

-14.1 

3 

1 .891 783E-12 

1 .89 1 783E- 1 2 

1.96 

0.06  ' 

8 

6.1 

8 

2.013388E-09 

2.0I3388E-09 

1.23 

0.00  ' 

9 

6.4 

8 

1.522937E-09 

1.522937E-09 

1.05 

0.32 

10 

6.8 

8 

6.1 18262E-10 

6. 1 1 8264E-09 

3.54 

0.78 

11 

7.1 

8 

1.663523E-10 

1.663523E-09 

3.95 

0.94 

12 

11.3 

8 

7.590633E- 1 1 

7.590636E-1 1 

4.01 

0.92  | 

13 

11.3 

8 

1.883713E-1 1 

1 .88371 3  E- 1 1 

4.16 

0.58  ! 

14 

11.3 

8 

3.359880E-12 

3.359880E-12 

4.35 

0.15 

15 

12.2 

8 

5.32301 8E-1 3 

5.32302 1 E- 1 3 

4.51 

0.90 

16 

12.2 

8 

6.841 38 1 E- 1 4 

6.84 1 383E- 1 4 

4.70 

0.65 

17 

13.1 

8 

8.156695E-15 

8. 156692E-15 

4.85 

0.89 

18 

25.6 

18 

5.861817E-15 

5.86 1 8 1 9E- 1 5 

4.61 

0.50 

19 

28.7 

19 

2.75041 5E-15 

2.75041 7E- 15 

4.42 

1.16 

20 

32.9 

20 

1 .38 1 999E- 1 5 

1.382000E-15 

4.22 

0.73 

21 

33.4 

20 

3.44841 9E- 1 6 

3.448422E-16 

4.07 

1.17 

22 

33.4 

20 

5.740436E-17 

5.740418E-1 7 

4.21 

1.71 

23 

34.5 

20 

8.338935E-18 

8.339004E-18 

4.72 

2.14 

Fig.  3.6.3:  Example  of  the  hybrid  algorithm  with  clustering  for  t  -  I  . 


example:  The  modified  hybrid  algorithm,  with  clustering,  is  illustrated  in  Fig.  3.6.3  for  a  col¬ 
lection  of  abscissae  which  includes  confluent,  close  and  well-separated  data  points. 
The  clustering  function  is  g(k)  =  k.  The  third  column  of  Fig.  3.6.3  indicates  the 
resulting  clustering  of  the  abscissae.  The  fourth  column  contains  the  top  row  of  the 
divided  difference  table  computed  in  single  precision  with  about  seven  decimal 
digits.  The  fifth  has,  for  comparison,  the  same  differences  computed  in  double  pre¬ 
cision.  Finally,  a  priori  error  bounds,  calculated  from  (3.4.15)  and  a  growth  factor 
bound  from  (3.2.8),  and  the  actual  relative  error  are  given  in  a  digits  lost  form. 
Complete  tables  corresponding  to  Fig.  3.6.3  are  presented  in  Appendix  D.  Fig.  3.6  4 
repeats  the  same  computation,  but  with  r-2.  Finally,  Fig.  3.6.5  shows  the  result  of 
computing  the  entire  table  in  one  scaling  and  squaring  for  A0nexp:. 


rAD-AQ8fl  751  CALIFORNIA  UNIV  BERKELEY  ELECTRONICS  RESEARCH  LAB  F/6  li 

ACCURATE  COMPUTATION  OF  DIVIDED  DIFFERENCES. <U> 

MAY  80  A  C  MCCURDY  N00014-76-C-0013 
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Divided  Differences  by  the  Hybrid  Algorithm 

n 

Abscissae 

f. 

Clustering 

Ao'expj  by 
algorithm 

Correct  values 
to  7  digits 

A  priori 
bound 

Relative 

error 

0 

-34.5 

0 

1.080639E-30 

1.080639E-30 

0.00 

0.00 

1 

-33.1 

1 

1.192152E-29 

1 .192 1 52E-29 

0.50 

0.00 

2 

-32.9 

1 

1 .9861 83E-29 

1 .9861 83E-29 

0.77 

0.26 

3 

-14.4 

3 

4.467966E-17 

4.467966E- 1 7 

0.51 

0.16 

4 

-14.4 

3 

8.233205E-I7 

8.233205E-17 

0.58 

0.00 

5 

-14.4 

3 

7.604 186E- 17 

7.6041 85E- 1 7 

1.85 

0.13 

6 

-14.4 

3 

4.692805E-I7 

4.692805E-17 

1.97 

0.08 

7 

-14.1 

3 

2.444065E-17 

2.444065E-I7 

2.03 

0.07 

8 

6.1 

8 

8.977370E-07 

8.977370E-07 

0.95 

0.00 

9 

6.4 

8 

1.963  505E-06 

1.963505E-06 

0.72 

0.00 

10 

6.8 

8 

2.309258E-06 

2.309258E-06 

2.38 

0.50 

11 

7.1 

8 

1.760567E-06 

1.760566E-06 

2.59 

0.32 

12 

11.3 

12 

1.305495E-05 

1.305495E-05 

2.23 

0.00 

13 

11.3 

12 

1.217449E-05 

1.217449E-05 

1.99 

0.53 

14 

11.3 

12 

6.558480E-06 

6.558482E-06 

3.41 

0.80 

15 

12.2 

12 

3.453429E-06 

3.453430E-06 

3.86 

0.43 

16 

12.2 

12 

1.281568E-06 

1.281569E-06 

4.26 

0.98 

17 

13.1 

12 

4.751205E-07 

4.751 204E-07 

4.48 

0.55 

18 

25.6 

18 

9.602053E-04 

9.602055E-04 

3.93 

0.63 

19 

28.7 

19 

1.412638E-02 

1.412640E-02 

3.28 

1.17 

20 

32.9 

20 

4.335487E-01 

4.335489E-01 

2.54 

0.93 

21 

33.4 

20 

6.106924E-01 

6.106929E-01 

1.92 

1.04 

22 

33.4 

20 

3.836830E-01 

3.836832E-01 

2.55 

0.84 

23 

34.5 

24 

2.38I056E-01 

2.38I055E-01 

2.94 

0.54 

Fig.  3.6.4:  Example  of  the  hybrid  algorithm  with  clustering  for  r  —  2 . 


Special  methods  for  low  order  differences.  It  is  sometimes  possible  to  compute  low  order 
divided  differences  by  a  special  formula.  From  (3.1.2)  where 


AiJexp, 


r<f  i+*n>/2  sinhfr  (g ,  -  £p)/2l 
*  (fi-fo>/2 


we  see  that  first  order  differences  may  always  be  computed  accurately  when  a  good  sinh  func¬ 
tion  is  available.  Error  growth  in  using  the  standard  divided  difference  formula  is  primarily 
dependent  on  errors  propagated  from  low  order  differences.  Special  computation  of  these 
differences  may  be  very  effective*  in  reducing  errors  in  higher  differences  and  in  extending  the 

tFig.  3.8.3  gives  an  example  (for  complex  abscissae)  of  dramatic  improvement  in  the 
error  when  first  order  divided  differences  are  computed  by  a  special  formula. 


area  over  which  this  simple  formula  may  be  used.  In  addition,  scaling  and  squaring  never  need 
be  used  for  these  low  order  differences. 


Divided  Differences  by  Scaling  and  Squaring 

n 

Abscissae 

Clustering 

Atfexpj  by 
algorithm 

Correct  values 
to  7  digits 

A  priori 
bound 

Relative 

error 

0 

-34.5 

0 

1.080639E-30 

1.080639E-30 

0.00 

0.00 

1 

-33.1 

0 

1.192152E-29 

1 . 1 92 1 52E-29 

0.50 

0.00 

2 

-32.9 

0 

1.986174E-29 

1 .986 1 83E-29 

4.85 

1.88 

3 

-14.4 

0 

4.467955E-17 

4.467966E-17 

4.85 

1.60 

4 

-14.4 

o 

8.233 186E- 17 

8.233205E-I7 

4.85 

1.57 

5 

-14.4 

0 

7.604169E-17 

7.6041 85E- 17 

4.85 

1.55 

6 

-14.4 

0 

4.692796E-17 

4.692805E-17 

4.85 

1.50 

7 

-14.1 

0 

2.444060E-17 

2.444065E-I7 

4.85 

1.56 

8 

6.1 

0 

8.977336E-07 

8.977370E-07 

4.85 

1.80 

9 

6.4 

0 

1.963496E-06 

1.963505  E-06 

4.85 

1.89 

10 

6.8 

0 

2.309243E-06 

2.309258E-06 

4.85 

2.06 

11 

7.1 

0 

1.760551  E-06 

1.760566E-06 

4.85 

2.16 

12 

11.3 

0 

1.305484E-05 

1. 305495 E-05 

4.85 

2.13 

13 

11.3 

0 

1.217439E-05 

I.217449E-05 

4.85 

2.12 

14 

11.3 

0 

6.558430E-06 

6.558482E-06 

4.85 

2.12 

15 

12.2 

0 

3.453403E-06 

3.453430E-06 

4.85 

2.11 

16 

12.2 

0 

1.28 1 559E-06 

1.281569E-06 

4.85 

2.11 

17 

13.1 

0 

4.75 1 1 75E-07 

4.751204E-07 

4.85 

2.02 

18 

25.6 

0 

9.601986E-04 

9.602055E-04 

4.85 

2.08 

19 

28.7 

0 

1.412635E-02 

1.412640E-02 

4.85 

1.71 

20 

32.9 

0 

4.335491  E-01 

4.335489E-01 

4.85 

1.00 

21 

33.4 

0 

6.106924E-01 

6.106929E-01 

4.85 

1.04 

22 

33.4 

0 

3.836829E-01 

3.836832E-01 

4.85 

1.13 

23 

34.5 

0 

2.381053E-01 

2.381055E-01 

4.85 

1.09 

Fig.  3.6.5:  Example  of  the  scaling  and  squaring  algorthm  for  r  —  2 . 


Second  order  differences  also  may  be  computed  accurately  by  a  special  formula  when  a 
routine  is  available  to  evaluate  the  function 


h(()  = 


e(—  1 


ao 

I 


O'+D! 


accurately  for  all  (.  Let  (0  <  f i  ^  then 
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1 

i 

I 


! 

i 

4 


—  reT(,[( 


r(f2-f|) 


1>  -  ( 


T«0-<|) 


-1)} 


Ter('[h[T(ti-(,)]  -  6[r(f0-f|))} . 


Because  sign{/>(€)] -sign({)  and  while  fo-fi<0,  the  subtraction  is  actually  an 

addition  of  non-negative  numbers. 

Perturbations  and  shifts  in  the  abscissae.  Abscissae  used  in  computing  divided  differences 
may  be  obtained  either  experimentally,  or  as  the  result  of  earlier  computations.  In  either  case 
we  may  be  uncertain  what  are  the  exact  abscissae  (represented  here  by  the  vector  x,  say).  The 
abscissae,  say  x,  we  have  in  hand  are  only  approximations.  The  most  we  can  expect  is  to  have 
a  bound  in  terms  of  x  on  our  uncertainty  in  the  value  of  x.  Thus  given  x  and  a  bound  on  the 
uncertainty,  we  ask  how  far  can  the  divided  difference  A"expr(x)  be  from  A"expr(x).  That  is, 
how  unsure  are  we  of  the  value  of  a  divided  difference,  given  our  doubt  about  its  data. 

As  an  example,  we  have  presented  without  comment  several  formulas  in  which  abscissae 
are  shifted  by  a  constant  amount,  say  a.  In  finite  precision  arithmetic,  a  computed  shifted 
abscissa  +  satisfies 

[/?(£  +  <*)  —  (f  +  a)|  ^  (|£|  +  |a|)«. 

To  have  a  uniform  bound  for  all  abscissae  represented  in  a  vector  x,  we  write 
|l//(x+cr«)  -  (x  +  a«)|U  <  (||x||„+  |a|)«  .f 

The  bound  describes  our  maximum  uncertainty  in  where  the  exact  shifted  vector  of  abscissae 
lies,  given  knowledge  only  of  the  computed  vector. 

It  is  convenient  to  think  of  x  as  a  perturbation  of  the  given  vector  x.  The  following  per¬ 
turbation  bound  describes  the  sensitivity  of  A  "exp,  to  a  bounded  change  in  its  abscissae. 


tRecall  that  u  is  a  vector  of  l’s,  u  •(!.  1 . 1). 
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Perturbation 

bound.  Let  Jr-  . . . 

be  a  vector  of 

abscissae,  and 

1 

© 

,f„)  a  perturbation  of  x  such  that 

max  <  y« 

for  a  constant  y. 

o 

Then 

|A"expr(i)  -  A"expr0r)|  <  ( er 

l)  A"expr(jf) . 

(3.6.2) 

proof:  From  Theorem  1  in  §3.2,  A"expr  is  increasing  in  each  of  its  abscissae,  thus 
A"expr(jr -yeu)  <  A"exp,W)  <  k"cxpT(x  +  yeu) . 

By  the  translation  property  (3.1.1), 

<?-n',.A',expTOr)  <  A"expr(x)  <  eri’‘  A"expr(x) .  □ 

For  small  rye  the  bound  (3.6.2)  is  equivalent  to  a  relative  error  of  size  rye.  Hence  com¬ 
putational  errors  may  be  viewed  in  the  same  way  as  uncertainties  in  the  data.  In  particular 
when  data  uncertainties  of  size  ye  lead  to  uncertainties  of  size  rye  in  the  value  of  A"expr(.v) 
relative  to  A"expT0r),  computational  errors  of  comparable,  or  smaller,  size  do  not  greatly 
increase  our  uncertainty.  Thus,  there  may  be  no  reason  to  compute  A"expr(jr)  to  greater  accu¬ 
racy  than  about  rye.  Hence  our  uncertainty  in  the  data  helps  answer  the  question  of  how  much 
accuracy  we  are  justified  in  demanding  when  computing  divided  differences.  We  may,  then, 
use  the  fast  standard  scheme  more  in  practice,  as  the  data  may  not  warrant  using  more  accu¬ 
rate,  but  more  costly,  methods. 

Additional  modifications  of  the  basic  hybrid  algorithm  may  be  desireable  in  practice.  For 
example  Ao'expT  decreases  as  r"/n\ ;  so  special  provisions  may  be  required  to  represent  small 
numbers  during  computation.  These  details,  however,  must  not  obscure  the  important  fact 
about  the  hybrid  algorithm,  which  is  that  real  exponential  divided  differences  can  be  computed 
with  high  relative  accuracy.  Such  a  general  statement  cannot  be  made  when  the  abscissae  are 
complex.  However,  a  hybrid  type  algorithm  with  error  bounds  comparable  to  the  above  can  be 
developed  for  some  arrangements  of  complex  abscissae.  We  turn  to  such  a  problem  in  the  next 
three  sections. 


§3.7 


9! 


3.7  Divided  differences  of  the  exponential  function  with  complex  abscissae. 

For  applications  exponential  divided  differences  with  complex  abscissae  are  more  impor¬ 
tant  than  the  real  case.  In  particular  a  real  non-Hermitian  matrix  A  may  have  complex  eigen¬ 
values.  These  eigenvalues  are  the  abscissae  used  in  forming  coefficients  of  the  Newton  polyno¬ 
mial  form  of  expire ).  Therefore  it  is  important  to  understand  which  aspects  of  the  real  case 
go  over  to  the  complex  case,  and  which  do  not. 

The  algorithms  presented  earlier  are  applicable  to  complex  abscissae.  The  theory  used  to 
derive  the  Taylor  algorithm,  scaling  and  squaring,  and  even  the  hybrid  algorithm,  depends  only 
upon  the  exponential  function  itself.  There  is  no  need  to  distinguish  between  real  and  complex 
data  points.  Our  error  bounds  and  decision  criteria,  however,  do  depend  explicitly  upon  the 
fact  that  real  exponential  divided  differences  are  positive.  Since  complex  differences  can  be 
zero,  we  must  abandon  the  idea  of  strict  relative  error  bounds.  Instead,  we  give  error  bounds 
relative  to  a  quantity  that  bounds  or  estimates  our  divided  difference. 

In  this  section  we  examine  a  few  special  cases  of  complex  exponential  divided  differences 
in  order  to  gain  a  better  understanding  of  the  behavior  of  such  differences.  In  particular  we 
shall  observe  how  these  divided  differences  are  affected  by  the  imaginary  parts  of  the  abscissae. 
Later  we  indicate  how  our  algorithms  may  be  applied. 

We  continue  studying  divided  differences  of  the  function  ,/  -expT  with  parameter  t  >0. 

Our  sequence  of  abscissae  Z  =  ({o.  Ci . {„,... )  may  now  contain  complex  elements.  We 

look  at  three  special  arrangements  of  the  abscissae.  (I)  The  abscissae  lie  on  a  line  in  the  com¬ 
plex  plane  and  are  evenly  spaced  along  this  line.  (2)  The  sequence  of  abscissae  consists  of 
repetitions  of  two  points,  {  and  -£;  we  also  look  at  the  case  where  the  two  points  are  conju¬ 
gates  £  and  £.  (3)  Finally,  we  examine  the  case  where  the  sequence  of  data  points  consists 
exclusively  of  conjugate  pair  points.  In  the  Rrst  two  examples  we  achieve  explicit  formulas  for 
the  divided  differences.  In  the  Anal  case,  we  characterize  the  differences  by  upper  bounds  on 
them.  This  final  case  is  of  most  interest  in  matrix  function  computations  because  the  eigen¬ 
values  of  real  matrices  are  either  real  or  members  of  complex  conjugate  pairs. 

Evenly  spaced,  linear  abscissae.  On  a  line  the  abscissae  can  be  ordered.  Let  (0  be  an  extreme 
data  point  and  let  25  be  the  spacing  between  the  abscissae.  Then 

Z  “{{<>•  {o+25. £o+4$ . £0+2n5,...)  is  the  sequence  of  data  points.  Exactly  as  in  the  real 

case  in  §3.2, 


(3.7.1) 


A,i'expr 


«! 


<>*-1 

28 


)" 


1  ^  sinh(r8) 

ill  8 


We  note  that  Id'exp,— 0  if,  and  only  if.  t8*  kni  for  some  non-zero  integer  k.  Thus  high  order 
divided  differences  can  be  zero.  Since  two  points  lie  on  a  line,  this  also  implies  that  first 
divided  differences  are  zero  if.  and  only  if,  their  abscissae  are  separated  by  k  iri. 

Suppose  8  is  pure  imaginary,  say  8  —  vi.  Let  us  observe  how  |Anexpr|  varies  with  v.  We 

have 


Ud'expT| 


sinh(ri/?) 

vi 


I  rtftj  sin(rv) 
n\  v 


where  |o=Re((n)-  |Ad'expT|,  then,  behaves  as  a  damped  sine  wave,  becoming  smaller  with 
increasing  v.  It  has  local  maxima  when  tan(T»») »rv.  For  fo**0  and  r-1,  the  table  in  Fig. 
3.7.1  lists  some  of  these  maxima  for  n  -  1.2 . 7. 


V 

|A<jexp| 

U<?exp| 

|A(Jexp| 

lioexpl 

|A<jexp| 

Ufcxp| 

|l07expl 

0 

1.00 

5.00E-1 

1.67E-1 

4.17E-2 

8.33E-3 

1.39E-3 

1.98E-4 

4.49 

2.17E-1 

2.36E-2 

1.71E-3 

9.28E-5 

4.03E-6 

1.46E-7 

4.53E-9 

7.73 

1.28E-1 

8.24E-3 

3.53E-4 

1.13E-5 

2.91E-7 

6.22E-9 

1.14E-10 

10.90 

9.13E-2 

4.17E-3 

1.27E-4 

2.90E-6 

5.29E-8 

8.06E-10 

1.05E-1 1 

14.07 

7.09E-2 

2.51E-3 

5.94E-5 

1.05E-6 

1.49E-8 

1.77E-10 

1.79E-12 

17.22 

5.80E-2 

1.68E-3 

3.25E-5 

4.71E-7 

5.46E-9 

5.27E-U 

4.37E-13 

20.37 

4.90E-2 

I.20E-3 

1.96E-5 

2.41E-7 

2.36E-9 

1.93E-1 1 

1.35E-13 

23.52 

4.25E-2 

9.02E-4 

1.28E-5 

1.36E-7 

I.I5E-9 

8.16E-I2 

4.95E-I4 

26.67 

3.75E-2 

7.02E-4 

8.77E-6 

8.22E-8 

6.I6E-I0 

3.85E-12 

2.06E-14 

29.81 

3.35E-2 

5.62E-4 

6.28E-6 

5.26E-8 

3.53E-IO 

1.97E-12 

9.44E-15 

Fig.  3.7.1:  Maxima  of  lAtfexpl,  as  a  function  of  v,  for  evenly 
spaced  imaginary  abscissae. 


The  magnitude  of  these  divided  differences  is  strongly  affected  by  the  difference  in  the 
imaginary  parts  of  the  abscissae.  Our  study  of  complex  exponential  divided  differences  must 
take  this  into  account.  The  next  example  even  more  clearly  illustrates  this  dependence  on  the 
imaginary  parts. 
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Two-point  exponential  divide#  differences,  in  an  example  in  §2.8  we  saw  that  divided 
differences  of  exponential  functions  for  sequences  of  data  points  like  Z I. 
where  a  point  and  its  negative  are  repeated,  have  many  special  properties.  In  particular,  we 
found  that  the  related  functions 


b„(0  =  -j{(  Jr  "exp, )({,-{ . 0  +  (AJ*exprX-s.£. 


a,M)  =  (42"+lexpr)({.-{ . £.-£) 


(3.7.2a) 


(3.7.2b) 


satisfy  the  recurrences 


bj()  - 


rb, ,_,(£)  -  (2n— l)fl„_|({) 
°M - - 


(3.7.3a) 


(3.7.3b) 


for  n  -  1 , 2 . where 


b0(0  “  cosh(r{) 


sinh(T{) 


(3.7.4a) 


(3.7.4b) 


From  these  relations,  we  show  that  the  functions  b„  and  a„  are  representable  in  terms  of  spher¬ 
ical  Bessel  functions,  commonly  denoted  j„.  In  addition,  we  derive  a  simple  assymptotic 
expression  for  the  two-point  divided  difference  (42"expr)(£,-{ . {)  as  r|{|— *°°. 


Representation  of  two-point  exponential  divided  differences.  For  each  n  -0, 1.2 . 


b"(°  m  2"nWO"-'J"~'UT° 

-»+ 1 

°M>  -  ■ 

where  the  j„  are  spherical  Bessel  functions.  Also  as  t|{|  —  «>, 
(AJ"expr)({,-{ . {)  -  -p^r{. 


(3.7.5a) 


(3.7.5b) 


(3.7.6) 


proof:  Spherical  Bessel  functions1  are  related  to  the  more  familiar  Bessel  functions  of  the  first 

tThe  introduction  to  the  National  Bureau  of  Standards'  Tables  of  Spherical  Bessel  Func¬ 
tions  (1947)  gives  a  brief  explanation  of  these  functions. 
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kind.  Jm%  according  to 

j„(u)  *■  yJir/2u>  J,l+i/.(u) . 

where  the  index  /w  — w+'/j  indicates  a  half-order  Bessel  function.  The  well-known  Bessel 
recurrence 


2nUlu(at)  ”  wf4f|lwl  +  /„_|(a>)) 


becomes 


(2n+l)j„(a>)  -  &»{>„+|(a»)  +>„_|(<o)) 


for  spherical  Bessel  functions.  Initially 


jo(u)  m  s‘n^M'  and  y_i (o*) 

<o 


cos(cu) 

it t 


(3.7.7) 


(3.7.8) 


When  fl—0  in  (3.7.5a-b),  a  comparison  of  (3.7.4a-b)  with  (3.7.8)  shows  that  initially  the 
j- 1  and  j0  in  (3.7.5a-b)  are  spherical  Bessel  functions.  For  general  n  we  derive  the  Bessel 
recurrence  (3.7.7)  for  j„  from  the  recurrences  (3.7.3a-b)  for  b„  and  a„.  Inserting  (3.7.5a-b) 
into  (3.7.3b)  yields 


T«+l  J 

2 "it!(/{)"-/"(,TC)  "  2*0  1  2"-,(/»— 1)!(/{)"-2 J'"“2(,ri> 


2"~' 


(2n—\)r" 

(«— 1)!(/()"-1 


y„_,(tr{)| 


or 

j„(irO  •  -  ,(2^-  -■->„.,(/>{) . 

When  this  is  rearranged  and  the  index  n  is  increased  by  1,  we  obtain  (3.7.7)  with  «  —  /t{; 
hence  each  j„  in  (3.7.Sa-b)  is  a  spherical  Bessel  function. 

For  large  |<u|  spherical  Bessel  functions  behave,  assymptotically,  as 

j„(u)  —  — cos[«- (/i+Dw/21 . 
u 

Thus  as  r|{j  —  oo  we  have 

<AJ"expr)({.-{ . 0  -  MO  +  CMO 

T#  +  l 

-  2votFrU-,(/TC)  ■  U'Ut0) 
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4 


r" 

—  j-^\cos(iTl-  iiit/2)  -  isin(/r£  —  «ir/2>) 


-£>*  □ 


2"«!{" 


The  translation  property  (3.1.1)  provides  an  immediate  corollary  to  the  above  when  the 
sequence  of  abscissae  consists  of  {  and  repeated. 


Corollary:  Let  the  sequence  of  abscissae  Z «{£.{. where  £-£  +  /•»)  and  its  conju¬ 
gate  {  «  i  -  iy  are  repeated.  Then  for  each  /;  -  0. 1 , 2 . 


,  \#l  —  l_ff+l 

Regexp,)  -  e"bjiv)  -  LJ±r_ertj  (_T  , 

2"n\y‘  1 

(3.7.9a) 

Im(A(J"expr)  -  -»j-AoJ"+lexpr  -  yeTfa„(iy)  -  ■~;I-TcT(j„(-ry) . 

2  n\y 

(3.7.9b) 

Further  as  tj  — <», 

4<?"expr  -  T--p^. 

2"niy" 

(3.7.10) 

proof:  By  the  translation  property. 

A0J"exPr  =  (42"expr)({,{,  ...,£)-  eT(  (^hlexp.)Uy,-iy.  .  . 

,iy)  • 

The  results  follow  from  (3.7.5a-b)  and  (3.7.6)  by  inserting  iy  for  £,  and  then  multiplying  by 

e7(.  □ 

From  (3.7.10), 

|4o"expr|  ~  ner(. 

2  n'.y 

(3.7.11) 

The  imaginary  part  leads  to  a  t damping  of  A(j''expT.  Also  since  (— t-tj)  and  >„(-■ r y)  are 

never  simultaneously  zero,  for  all  r  >  0  the  divided  difference  4<?"expT  ^  0. 

Exponential  divided  differences  with  conjugate  pair  abscissae.  We  now  turn  to  the  case  of  a 
sequence  of  abscissae  consisting  of  conjugate  pair  elements.  In  particular  let 
Z*{Co.Co.(i.Ci . I  where  and  with  each  y,  >0.  The  fol¬ 

lowing  bounds,  depending  on  divided  differences  of  expr  for  the  real  abscissae  f ,  alone,  help  to 
describe  the  dependence  of  conjugate  pair  exponential  divided  differences  on  both  the  real  and 
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imaginary  parts  of  the  abscissae. 


Bounds  for  conjugate  pair  exponential  divided  differences.  Let  Z  - 1£„.  ■ 

be  a  sequence  of  conjugate  pair  abscissae.  Then  for  each  n  >  0, 

|AoMexpr|  ^  (fi^-'^exp. 

(2.7  12a) 

and 

Uo"+lexp,|  <  <riny)"‘  AfflexpT. 

(2.7.1 2b) 

/- 0 

proof:  The  proof  is  by  induction  on  />.  We  note  first,  employing  a  remark  after  (2.1.3),  that 
since 

,  (A2"exp,)({0.?o.  •  •  •  ~  U2"exPr)({0.?o . C„_,. 

An  exp,  - - - - = - 


-  — Im(l^"expr) , 

Vn 


(3.7.12b)  is  an  immediate  consequence  of  (3.7.12a).  When  n  -0  and  r  >  0, 

|A0°exprl  -  kTt°l  -erfn-l(°0exPr. 

and  (3.7.12a)  certainly  holds.  Now  let  us  assume  it  is  true  for  all  orders  up  to  (2/r — 2).  That 
is,  we  assume  for  all  r  >  0 


|A(?"~2expr|  <  (rjfr/r'  Ar'exp,, 
;-0 


and  hence 


«-i 


Uo"_lexPrl  < 

7*0 


By  the  recursive  integral  formula  (3.1.5), 

T 

Atf"exp,  -  er<"J* . 


Thus 
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|A,t"expr|  ^  eH“ JV,r<'iA,?"~lexp(r|  da 
n 

M  - 1  f 

<  (nv.>"  rT("f  <'  ,i"'Af"'lexp,r  da 
i-o  n 

-  (nn/)',  A/oexpr.  □ 

/» 0 

When  the  real  parts  £,  of  the  abscissae  are  equal,  that  is  • 

•  -  f„,  this  corollary 

follows. 

Corollary:  Let  Z  «■» |f  ±  /tj,,  j  *0,  .  .  .  .n ).  Then  for  each  n  >  0. 

|Ao"expr|  <£ 

i-0  n  • 

(3.7.13a) 

and 

!lo"+,expr|  < 

1-0  "• 

(3.7.13b) 

With  the  exception  of  the  factor  2",  the  bounds  in  our  corollary  resemble  our  assymptotic 
results  for  two-point  conjugate  pair  divided  differences.  This  leads  us  to  suspect  that 

2-"(nn/>-'-Kexpr  and  2“"(]X>jy)'l-A/)1expr  (3.7.14) 

i-O  j- 0 

are  reasonable  estimates  of  !Ao"expT|  and  |A<j"+,exp.|,  respectively,  when  the  -rj,  are  large. 
The  values  in  Fig.  3.7.2  illustrate  this.  Note  that  not  every  estimated  value  by  (3.7.14)  is  large 
enough  to  be  a  bound. 

General  complex  exponential  divided  differences.  When  we  are  unable  to  make  assumptions 
about  the  abscissae,  we  can  say  little  about  the  behavior  of  the  divided  differences.  Just  as  the 
simple  bound  derived  from  (2.1.12), 


J Ao'expr|  <  max  \r"e 


\ 

i 


poorly  describes  the  behavior  of  4<j'expr  when  is  large,  even  when  all  the  abscissae  are 

real,  a  bound  depending  only  on  the  real  parts  of  complex  abscissae  poorly  describes  |Ao'expr| 
when  some  data  points  have  large  imaginary  parts.  All  complex  exponential  divided  differences 
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n 

c. 

Ud'expl 

bounds  from 
(3.7.12a-b) 

estimates  from 
(3.7.14) 

n 

(0.+  10/) 

1.00 

1.00 

1.00 

i 

(0,-10/) 

5.44E-2 

1.00E-1 

I.00E-1 

2 

(l,+9/> 

8.39E-2 

1.72E-1 

8.59E-2 

3 

(1,-9/) 

9.31E-3 

1.9IE-2 

9.55E-3 

4 

(2. +  11/) 

3.64E-3 

1.64E-2 

4.10E-3 

5 

(2,-11/) 

2.71E-4 

1.49E-3 

3.73E-4 

6 

(3, +  10/) 

9.18E-5 

8.54E-4 

1.07E-4 

7 

(3.-10/) 

3.16E-6 

8.54E-5 

1.07E-5 

8 

(4, +9/) 

2.41E-6 

3.67E-5 

2.29E-6 

D 

(4,-9/) 

2.50E-7 

4.08  E-6 

2.55E-7 

|fl 

(5, +  10/) 

4.20E-8 

1.40E-6 

4.38E-8 

(5,-10/) 

1.96E-9 

1.40E-7 

_ 

4.38E-9 

Fig.  3.7.2:  Bounds  and  estimates  for  |Ad'exp|  with  conjugate  pair  abscissae. 


do  satisfy  the  following  bound  regardless  of  the  imaginary  parts  of  the  data  points. 

Upper  bound  on  |Aoexpr|  for  complex  abscissae.  Let  Z-{4o.£i . £„ . }  be  a  se¬ 
quence  of  complex  valued  abscissae,  and  let  £ ,  *•  Re(£;)  for  each  j  -  0. 1 . n .  Then 

|i<j'expT|  ^  VoexpT  (3.7.15) 

for  all  n  ^  0. 


proof:  Directly  from  (3.1.3),  namely 


f  1 

A<j'expT-J7 


o  o 


J*  exp [t£0  +  (£ i “Co) * i  +  '  •  •  d<r„  '  '  •  d<r2d<r 


we  have 


|Ao"expr|  <//•••  J  exp[rf0+(^|-fo)o,i+  •  '  '  +  (£„-f„-i )(rll]d<r„  ■  ■  •  </<r2</ a 

0  0  0 


-A^;expT.  a 


Comparing  (3.7.15)  with  our  conjugate  pair  bounds  shows  that  the  imaginary  parts  of  the 
abscissae  may  be  very  important  and  should  be  reflected  in  any  bounds  we  use.  In  the  next 
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sections  we  apply  our  divided  difference  algorithms  to  conjugate  pair  abscissae  and  use  the 
upper  bounds  and  estimates  we  have  presented  here  to  derive  error  bounds  for  the  computa¬ 
tions. 
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3.8  Divided  difference  tables  with  conjugate  pair  abscissae. 

The  elements  of  divided  difference  tables  whose  abscissae  are  conjugate  pairs  are  quite 
special  in  that  some  are  real  and  some  are  conjugates  of  others.  These  properties  do  not 
depend  upon  the  exponential  function,  but  apply  to  any  function  /  symmetric  about  the  real 
axis.  For  this  reason  we  digress  in  this  section  from  our  study  of  exp,  and  revert  to  the  func¬ 
tion  ./;  where  /({)—/({).  The  results  are  more  general  and  the  notation  is  more  compact. 
Applications  to  expr  follow  in  the  next  section. 


./■(C-j)  All/' 

Aij / 

aV 

/<{-*> 

A  -iJ 

A  lit 

/<{->) 

Ai|./ 

/({-o> 

A  V 

A  V 

AV 

A  V 

Aij. J 

A  V 

A V 

A  V 

aJi  j 

Aii/ 

AV 

AV 

Ala/' 

A  V 

AV 

AV 

/(Co) 

A  <)/ 

A<J/ 

Aj )j 

/(Ci) 

A| If 

A  if 

mi) 

A  2 If 

mi) 

Fig.  3.8.1:  Rearranged  divided  difference  table  A/ for  {{-j.C_2.C-1.t-o.Co.C1.C2.C3l- 


We  have  seen  that  abscissae  should  be  ordered  so  that  close  values  are  adjacent  to  each 
other.  It  follows  that  {  should  not  be  adjacent  to  {  when  lm({)  is  large,  as  would  be  natural.  A 
good,  but  unorthodox,  ordering  for  complex  conjugate  pairs  of  abscissae  is 

Z-{{, ,.{«-! . C1.C0.C0.C1 . {„}•  Some  extra  dividends  follow  from  this  choice  as  we 

show  below.  In  order  to  maintain  (as  closely  as  possible)  our  notation  A*/ to  indicate  the  use 

of  {,,{,+1 . {/+*.  we  write  {_,  for  {>•  The  table  in  Fig.  3.8.1  shows  a  typical  A/  where 

n  —  3.  The  entries  corresponding  to  the  top  row  in  a  naturally  ordered  table  are  underlined. 
These  are  the  entries  that  are  used,  for  example,  as  coefficients  in  a  Newton  polynomial. 

Let  Z  be  the  'step  matrix”  associated  with  the  sequence  Z,  that  is 


•  fc* 


i. 


I 


l 

¥ 

J 


/ 

n 
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C-.  1 


z  - 


C-n  1 

Co  1 

Ci 


I 

c„ 


This  extra  property  referred  to  above  is  that  both  Z  and  A/-/(Z)  are  Hermitian  about  the 
secondary  diagonal  (bottom  left  to  top  right).  This  property  is  most  conveniently  expressed  in 
terms  of  the  permutation  matrix 


Secondary  symmetry  of  If.  If  Z-(C« . {0.Co . CJ  and  /({)«/( £),  then  both 

/•Z and  /'A/are  Hermitian. 


proof:  Any  complex  matrix  5  is  Hermitian  if  IF  —  B.  We  denote  the  conjugate  transpose  of  /? 
by  fl*.  By  inspection,  l-Z  is  Hermitian, 

*  —  * 

/Z-Z*/. 

From  the  conjugate  transpose  property  (1.1.6), 

/(Z’)-/(Z)*; 
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so  employing  (1.1.3)  for  similarity  transformations. 

b/Ozrh 

-  h/O  iz)  -  /  /(Z). 

and  l-f(Z)  is  Hermitian.  □ 

This  result  means  every  divided  difference  lying  on  the  secondary  diagonal  of  A/ is  real, 
and  every  divided  difference  below  this  diagonal  has  a  conjugate  above  the  diagonal.  For  exam¬ 
ple  in  Fig.  3.8.1,  A lu/;  kl\f,  kl-J and  llj/are  all  real,  and  kljf  - kjf while  l-j/’-li 
Only  the  portion  of  the  table  on  and  below  the  secondary  diagonal  ever  need  be  computed.  For 
example,  and  all  differences  upon  which  it  depends  might  be  computed  by  a  series  method 
because  the  abscissae  may  be  close.  The  standard  formula  and  taking  conjugates  will  fill  out  the 
rest  of  the  table.  The  idea  is  illustrated  in  Fig.  3.9.1. 

From  our  discussion  here,  the  reordered  table  is  clearly  ideal  for  computation  by  a  hybrid 
algorithm.  We  consider  this  for  /-exp,  in  the  next  section. 

example:  Fig.  3.8.2  shows  that  reordering  abscissae  and  computing  first  order  differences  by  a 
special  formula  may  have  a  dramatic  effect  on  error  propagation  when  the  standard 
scheme  is  used.  The  abscissae  here  are 

—  —1-414214  ±  /8.585786 
{*,  -  1.412799  ±  dl.41563 
C±2-  1.414214  ±  dl.41421 
C*3  -  1.417039  ±  (11.41138 

First  order  (initial)  differences  were  computed  correct  to  seven  decimal  digits.  The 
standard  scheme  was  employed,  thereafter,  in  greater  precision  to  isolate  propagated 
errors.  The  figure  compares  divided  differences  trom  the  top  row  of  the  table  for 
the  natural  ordering  of  the  data  points  with  the  identical  differences  when  the  data 
points  are  reordered  as  suggested  in  this  section. 

Reordering  permits  many  differences,  for  which  error  growth  would  be  large  by  the  stan¬ 
dard  scheme,  to  be  computed  by  a  special  method.  We  see  here  with  reordering  that  close 
abscissae  contribute  only  to  first  and  second  order  differences.  These  first  order  differences  do 
not  contribute  to  error  growth  when  computed  specially.  However,  failure  to  compute  first 
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-t  r 


Order 

n 

Correct  values  to 

7  digits 

Natural  ordering  with 
special  computation  of 
first  differences 

0 

(-1.624537E-1,  1.808715E-1) 

(-1.624537E-1,  1.8087I5E-I) 

1 

(  2.I06638E-2,  0.0  ) 

(  2.106638E-2,  0.0  ) 

2 

(-5.269139E-2,  I.213605E-2) 

(-5.2691 40E-2,  1.213604E-2) 

3 

(  1.063 108E-3,  0.0  ) 

(  1.063107E-3,  3.965390E-I0) 

4 

(-1.1 14105E-4,  1.907577E-3) 

(-1.113 896E-4,  1.907435E-3) 

5 

(  1.671230E-4,  0.0  ) 

(  1 .67 1 291  E-4.-2.643648E-9) 

6 

(  3.230809E-S,  1.838758E-S) 

(  3.215634E-5,  1.820492E-5) 

7 

(  1.61 1337E-6,  0.0  ) 

(  1 .547398E-6.-3.59791 4E-8) 

Order 

n 

Reordering  with 
special  computation  of 
first  differences 

Reordering  without 
special  computation  of 
first  differences 

0 

(-1.624537E-1,  1.80871SE-1) 

(-1.624537E-1,  1.808715E-1) 

1 

(  2.106638E-2.0.0  ) 

(  2.106639E-2,  0.0  ) 

2 

(-5.269141  E-2,  1.213605E-2) 

(-5.269139E-2,  1.213604E-2) 

3 

(  1.063108E-3,  0.0  ) 

(  1.063108E-3.0.0  ) 

4 

(-1.1141 07E-4,  1.907577E-3) 

(-1.11 3844E-4,  1.907667E-3) 

5 

(  1.671 230E-4,  0.0  ) 

(  1.671309E-4,  0.0  ) 

6 

(  3.230604E-5,  1.838538E-5) 

(  3.203465E-5,  1. 80291 5E-5) 

7 

(  1.61  U44E-6,  0.0  ) 

(  1.579928E-6,  0.0  ) 

j 


Fig.  3.8.2:  Effects  of  reordering  data  points  and  special  computation  of 
first  divided  differences  on  A  "exp. 

order  differences  accurately  destroys  any  benefit  from  reordering,  as  the  numbers  show. 


When  clusters  of  close  abscissae  are  small,  as  here,  reordering  the  abscissae  makes  special 
computation  of  low  order  differences  very  effective  in  controlling  error  growth.  In  the  next  sec¬ 
tion  we  shall  see  that  a  hybrid  algorithm  effects  even  more  dramatic  improvements  in  accuracy. 
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3.9  A  hybrid  algorithm  for  Aexp,  with  conjugate  pair  abscissae. 

The  table  rearrangement  presented  in  the  last  section  strongly  suggests  implementing  a 
hybrid  algorithm  for  computing  Aexp,  when  the  abscissae  are  conjugate  pairs.  A  hybrid  algo¬ 
rithm  using  scaling  and  squaring,  as  well  as  the  standard  scheme,  is  most  accurate  for  abscissae 
having  imaginary  parts  nearly  equal  in  absolute  value,  but  large.  Fig.  3.9. 1  illustrates  the  conju¬ 
gate  symmetry  relationships  in  a  reordered  conjugate  pair  divided  difference  table.  The  "s"  indi¬ 
cates  an  element  which  may  be  computed  by  scaling  and  squaring,  "x"  by  the  standard  scheme, 
while  "r"  means  the  element  is  real. 

s  x  x  x  r 

s  x  x  r  x 

s  x  r  x  x 

s  r  x  x  x 

s  s  s  s 

s  s  s 

s  s 
s 

Fig.  3.9.1:  Relation  of  entries  in  conjugate  pair  table. 


s  s  s 
s  s 
s 


AexpT  - 


13 

Aio 

15 

a5 

A4 

13 

a3o 

15 

a4 

Ai2 

13 

15 

A4 

Ai, 

Ai, 

*,{® 

44 

Aio 

Aio 

Aio 

/to 

Ao 

A<? 

A*1 

e*1' 

A,' 

A|J 

e** 

a2< 

er<* 

Fig.  3.9.2:  Conjugate  pair  divided  difference  table  showing  symmetries. 
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Algorithm:  Hybrid  algorithm  for  conjugate  pair  abscissae.  For  conjugate  pair  abscissae 

(,«£,  +  /'n,  anu  -/>/  with  tj,  »0,  y-0. 1 . /»,  form  the  divided  difference 

table  matrix  as  follows: 

1.  Reorder  the  sequence  of  data  points  as  {£- . 5-i.{-o.{o.(i . {J  and.  if  possible 

(by  reindexing  if  necessary),  so  that  f  i  <  ■  ■  •  <  f„. 

2.  Compute  Ao'exp,,  and  hence  each  A,*ex pT  for  /  —  0.1 . /;  and  A  -0. 1 . n-i,  by 

scaling  and  squaring  (§3.4). 

3.  For  each  /  -  0, 1 . n  compute  (-i  +  lsO  when  /  —  0 ) 

a.  on  the  secondary  diagonal  of  the  table 

Ai'+1expT  -  —  Im(Al(+iexpT) ; 

V, 


t 


b.  for  k  -  2/+2 . /+«+ 1  each  Af ,expr  by  the  standard  scheme,  e.g. 

Air+iexpT  -  Air 'ex  PT 


Ai,expT 


~  C-, 


4.  Fill  the  remainder  of  the  table  using  conjugate  symmetry  about  the  secondary  diagonal. 


When  n-3  the  matrix  in  Fig.  3.9.2  illustrates  the  relation  between  various  elements  in  a 
table.  (Some  references  to  the  function  expr  are  suppressed.)  In  the  hybrid  algorithm  entries 
below  the  horizontal  line  in  Fig.  3.9.2  are  computed  by  scaling  and  squaring  (step  2),  while 
entries  to  the  left  of  the  vertical  line  are  just  conjugates  of  these,  as  indicated. 


Next  in  step  3a  of  the  algorithm. 


Aioexpr 


expT({o)  ~  expr({0) 
(o-(o 


— Im[expr({<>)j 
Vo 


r(nsin(rv0) 

Vo 


This  and  the  row  immediately  below,  already  known  from  the  scaling  and  squaring  step,  permit 
completion  of  the  -0  table  row  by  the  standard  scheme  (step  3b).  Then  again  in  step  3a. 

Ai(expr  -  —  Im(Aioexpr) , 

V\ 

and  the  elements  to  the  right  of  Aiiexp,  are  computed  by  the  standard  formula  (step  3b).  For 
example 


Step  3  is,  then,  repeated  until  Al^exp,  is  computed.  The  remaining  elements  to  the  right  of 
the  vertical  line  are  conjugates  of  elements  above  the  horizontal  line,  as  indicated.  The  ele¬ 
ments  on  the  secondary  diagonal  (indicated  by  underlining  in  Fig.  3.9.2),  computed  in  step  3a, 
are  all  real. 


J_A  3 

_L_  A  3 

1  k3 

1  Ai 

no  0 

noni  fo 

V0VIV2  ffl 

nonin2n3  f" 

_!_  A  2 

1  A 2 

1  1} 

noAf° 

VoVi  <0 

V0V1V2  (n 

V0V1V2  f" 

JLaj 

_L -At 

A  2 

1  A  3 

no  (° 

noni  f° 

noni  fn 

noni  f" 

JL/*« 

no 

jla, 

noA<° 

no 

J-A4 

no  1 

K 

H 

K 

eT(' 

K 

er(> 

Fig.  3.9.3:  Table  of  upper  bounds  based  on  real  divided  differences. 


Upper  bounds.  For  an  error  analysis  of  the  hybrid  algorithm  we  must  first  develop  error  bounds 
on  the  scaling  and  squaring  portion  of  the  computation.  Then  we  can  see  how  these  errors  are 
propagated  during  the  remainder  of  the  computation  by  the  standard  formula.  An  examination 
of  the  upper  bounds  (3.7.12a-b)  and  (3.7.15)  yields  quantities  relative  to  which  we  may  con¬ 
struct  error  bounds.  For  example,  the  table  in  Fig.  3.9.2  is  bounded,  element  by  element,  by 
the  table  in  Fig  3.9.3.  Here  we  omit  the  function  reference  expr,  for  clarity,  and  point  out  that 

the  divided  differences  in  Fig.  3.9.3  are  for  the  real  abscissae  |fo.fi . f  »)•  Our  error 

bounds  will  be  relative  to  an  upper  bound  matrix  such  as  the  one  in  Fig.  3.9.3. 


A  ; 

J_4  J 

1  A  ' 

*  A  1 

1  A«  j 

Cll 

2r»n  f" 

*Vv,Vl 

8-nnm-n:  in 

&VoVlV2V> 

A  2 

1  A 2 

1 _ ^ : 

1  A 

2-no  f" 

4t,oT>.  f" 

4wiin:  e" 

ZVoV'.V: 

At 

J_A. 

1  A  1 

•  A  } 

1  A’  1 

2-no  f" 

2nom  f" 

4r,oT»!  f° 

4hom  \ 

Vo 

-La' 

2hoAf" 

J_a  ■’  1 
2-no 

er(n 

K 

L. 

Fig.  3.9.4:  Table  of  estimated  absolute  values. 

When  the  tj,  are  nearly  equal,  but  large,  the  error  may  be  estimated,  element  by  element, 
by  (k  times  the  matrix  in  Fig.  3.9.4.  The  constant  k  depends  only  on  errors  introduced  in  the 
scaling  and  squaring  part  of  the  algorithm.  This  result  is  not  surprising.  When  the  17 ,  are  iarge. 
the  standard  scheme  is  employed  only  for  well-separated  abscissae.  From  our  earlier  studies 
there  is  little  error  growth  in  this  case.  Indeed,  such  separation  of  the  data  points  is  the  reason 
for  reordering  them,  in  the  first  place. 

example:  Witn  the  data  from  the  example  at  the  end  of  §3.8.  namely 

i±0-  -1.414214  ±  /8.585786 
{±)  -  1.412799  ±  /1 1.41563 
-  1.414214  ±  / 1 1.41421 
C£3  -  1.417039  ±  (11.41138. 

the  tables  in  Fig.  3.9.5  show  that  upper  bounds,  as  in  Fig.  3.9.3,  and  estimated  abso¬ 
lute  values,  as  in  Fig.  3.9.4,  describe  the  size  of  the  divided  differences.  From  sym¬ 
metry,  only  the  portion  of  each  table  on  and  below  the  secondary  diagonal  is  shown 
The  divided  differences  themselves  are  listed  in  Fig.  3.9.6. 

Scaling  and  squaring  error  bounds.  Error  bounds  from  our  earlier  analysis  of  scaling  and 
squaring  in  §3.4  carry  over  immediately  to  that  portion  of  the  conjugate  pair  table  computed  by 
this  method.  The  bounds  are  no  longer  valid  relative  to  the  computed  difference  itself,  but 
rather  to  an  appropriate  upper  bound  on  this  difference.  A  quick  reexamination  of  the 


¥ 


13: 


I 


Layout  of  tables 

Correct  absolute  values 

• 

• 

Aijexp 

1.61E-6 

Aijexp 

A^xp 

• 

I.67E-4 

3.72E-5 

Ai,exp 

Aiiexp 

Aiiexp 

1.06E-3 

I.91E-3 

7.80E-4 

Aioexp 

Aioexp 

Aioexp 

Aioexp 

2.1  IE- 2  5.4IE-2 

4.26E-2 

1.74E-2 

Afjexp 

A<?exp 

Aoexp 

2.43E-I  1.08 

8.54E-1 

3.53E-1 

,<• 

A /exp 

A  t2exp 

4.11 

4.11 

2.06 

*** 

Aj'exp 

4.11 

4.12 

4.12 

Upper  bounds  (Fig.  3.9.3) 

Estimated  values  (Fig.  3.9.4) 

■ 

• 

3.01  E-5 

• 

> 

3.76E-6 

8.67E-4 

3.43E-4 

• 

2.17E-4 

4.29E-5 

■ 

1.39E-2 

9.90E-3 

3.92E-3 

6.97E-3 

2.47E-3 

9.80E-4 

2.83E-2 

1.59E-1 

1.13E-1 

4.47E-2 

2.83E-2  7.96E-2 

5.65E-2 

2.24E-2 

2.43E-1 

1.37 

9.70E-1 

3.84E-1 

2.43E-1  1.37 

9.70E-1 

3.84E-1 

4.11 

4.11 

2.06 

4.11 

4.11 

2.06 

4.11 

4.12 

4.11 

4.12 

4.12 

4.12 

Fig.  3.9.5:  Example  of  bounds  and  estimates  for  conjugate  pair 
divided  differences  in  Aexp . 


derivation  of  scaling  and  squaring  error  bounds  will  show  this.  We  study  only  the  double  preci 
sion  accumulation  case,  as  the  argument  is  exactly  the  same  in  the  single  precision  case. 
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Scaling  and  squaring  bound.  Consider  only  the  non-conjugate  portion  of  the  sequence  of 

abscissae,  namely  |{0.{i . {„).  For  double  precision  accumulation  of  inner  products, 

l/MAexpr)  -  AexpT|  <  *1*2*0 -  ll  Axexpr,  (3.9.1) 

where  *2  —  8.3259,  0  is  the  maximum  spread  in  the  abscissae,  and  Axexpr  is  the  related  di¬ 
vided  difference  table  whose  abscissae  are  X -  i  Re(£,),  j  —  0 . n }.  For  single  precision 

accumulation, 

L/?(lexpr)  -  AexpT|  <  «(«+l)[*|T0-  1]-Axexpr,  (3.9.2) 

where  *t  —  21.2950. 


proof:  We  first  compute  a  scaled  divided  difference  table  by  the  Taylor  algorithm.  The  expan¬ 
sion  point  a  may  be  the  center  of  the  smallest  circle  enclosing  the  data  points,  and  the  spread  0 
is  the  diameter  of  that  circle.  Let  the  data  points  be  ordered  so  that  the  real  parts  satisfy 
•••<*„.  Because  exponential  divided  differences  with  real  abscissae  are  increasing 
in  each  abscissa,  we  have 

r“eT(o 

~7T~  <  A4«P* 


and 


r"|gra|  <  T^_  r<f0+*/2> 

n\  "  n\e 


eTH/2'J~7i~  *  erW/2^foexpr. 


In  Appendix  B  we  derive  the  error  bound 

l/72Ud'expT)  -  Ao'expr|  <  t(2  +  T9/2)eT*tlT  ^  '  ■ 

n : 


Therefore, 


l^2(Aoexpr)  -  Id'exPrl  <  <(2  +  T0/2)er*  A{J)expT.  (3.9.3) 

The  Taylor  series  error  bound  (3.9.3)  applies  to  every  element  of  the  divided  difference 
table.  The  error  in  the  original  scaled  matrix  in  scaling  and  squaring,  then,  must  satisfy  the 
matrix  inequality 


L/?2(Aexpr,r)  -  Aexpj.J  <  «(2  +  2-(>+l,T0)er/r#Axexpr,r. 


(3.9.4) 
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where  2~'  is  the  scaling  factor. f  The  subscript  X  indicates  the  divided  difference  table  matrix 
Axexpr,r  has  real  abscissae  X  -  tin.  (i . (J- 

For 


/3  j  =  (2  +  2-<'+"r#)eJ',r*. 

the  error  in  each  element  of  /fj(Aex prjlr)  is  bounded  by  c/8,  times  the  corresponding  element 
of  Axexp2-)f.  This  ,  is  exactly  that  used  earlier  in  §3.4.  For  any  complex  matrix  B 

[ /f2(fi2)  -  B2\  <  c|«|2. 


Also  bound  (3.7. IS)  yields 


Uexpr,T|  <  lxexpr,r; 


hence, 

|Aexpr<JJ  <  Axexpj_().|iT . 

The  same  argument  that  led  up  to  (3.4.9)  gives 

L/?2(AexpT)  -  Aexpr|  <  «(2'(/8,  +  l)  -  l]AxexpT, 

where  2'($,  + 1)  —  1  is  the  same  as  in  (3.4.9).  It  is  minimized  in  the  same  way.  For  j  the  smal¬ 
lest  non-negative  integer  satisfying  (3.4.10),  namely 


2~'t9  <  1.3292, 


we  obtain 


L/f2(AexpT)  -  iexpT|  <  c[K2r0- l]  AxexpT 

where  «c2- 8.3259.  The  same  argument  shows  the  single  precision  bound  is  (3.4.15).  □ 


tThe  matrix  bound  (3.9.4)  does  not  hold,  rigorously,  when  Aexp2_,r  is  backfilled  from 
its  top  row.  This  is  because 

At", 'exp,  <  Ae‘-'exp,+  |{,+*-{,|-Af*exp, 

when  tj,+h  pa  rj,.  If  the  Taylor  formula  is  used  on  the  entire  table,  (3.9.4)  does  hold. 


j 
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Relations  (3.9.1)  and  (3.9.2)  mean  that  the  error  in  complex  exponential  divided 
differences,  relative  to  corresponding  real  divided  differences,  has  the  same  bound  for  scaling 
and  squaring  as  in  the  exclusively  real  case. 

Standard  scheme  error  bounds.  When  all  the  imaginary  parts  tj,,  j  -0, 1 . »,  are  large  and 

nearly  equal,  the  portion  of  the  table  computed  according  to  the  standard  scheme  satisfies  the 
following  error  bounds. 

Standard  scheme  error  bound.  When  each  t>,,  j  —0, 1 . n,  is  large  compared  with  /i/r, 

t/(Ai**+fexpr)  -  Ait+'expr|  <  «K(J][d/)''-Afn+/"lexpT  (3.9.5) 

/-  o 

for  each  k  —  0, 1,  ....  n  and  /*■  1,2 . n—k+ 1.  k  is  one  of  the  scaling  and  squaring  er¬ 

ror  coefficients  k2t9  —  1  or  (/7+1)[k|T0-  1],  depending  on  the  arithmetic  used. 


The  recursive  nature  of  the  standard  scheme  makes  it  easiest  to  describe  bounds  on  the 
propagation  and  growth  of  errors  in  terms  of  examples.  Also,  this  will  make  clear  what  large 
compared  with  n/r  means.  Errors  introduced  in  the  scaling  and  squaring  portion  of  the  compu¬ 
tation  of  the  table  in  Fig.  3.9.2  are  propagated  during  the  computation  of  the  remaining 
differences.  We  bound  these  errors  relative  to  the  table  in  Fig.  3.9.3. 

From  (3.9.1-2),  depending  on  our  arithmetic  assumptions, 

)j7(Aoexpr)  -  <MexpT|  <  e*c-A{‘0expr 

for  each  k  —0, 1 . n.  To  keep  the  analysis  simple  we  forget  that  all  zeroth,  and  even  first, 

order  differences  may  be  computed  specially  with  a  smaller  error  coefficient  than  k.  The 
difference  loexpT  =  <?  c°  in  Fig.  3.9.2  is  computed  with  an  absolute  error  5 §  =  fl{eTla)  -  eTln 
such  that  [fiol  ^  €Ker(°.  Now, 

Aloexpr  -  —  ImU^exp,); 

■90 


I  Sic 


<  «K  — 


r(0 


90 


the  propagated  error  8io  satisfies 


Since  A(jexpr  is  computed  by  scaling  and  squaring,  its  absolute  error  »d  satisfies 
j8o I  ^  ««  Af(|expr.  By  (2.4.3)  the  propagated  absolute  error  in  the  computed  Ai^exp,  is 


S 


2 

-0 


8<i  ~8-o 

Ci-C-o  ’ 


thus 


A  - 


|6in 


<  tK~ 


1  T(n 

««e*  Pr  +  f 

VO 


ki 


ol 


1)0+  1/t  k  ,  I 

T? — : — r« — AiexpT 
ICi-C-ol  ijo  ft 


where  bound  (3.2.3)  is  used  on  er(°.  When  i)o  =  i)i  and  Vo  >s  large  compared  with  1/t, 

i)o+  l/r  _  \ 
ki-C-ol  2  ’ 


or  is  even  smaller  than  '/:  when  the  difference  in  the  real  parts  f  |  -£_0  is  large.  6 £0  satisfies  the 
simple  bound 


|8i0|  =  L/fUioexp,)  -  Aiocxprl  <  e— Af‘flexp7. 

1)0 


One  more  step  makes  the  general  case  (3.9.5)  clear.  Since  [So (  ^  eK-A^exp.  by  scaling 
and  squaring. 


«K" 


A/nexpT  +  ““Af  exp, 
^0 _ 

l?2-4-ol 


^  1)0+ 2/t  k 

<  17 - r-T-« — A;oexp,. 

1C2-S-0!  no  0 


Thus  when  i)o  =  i)i  =1);  and  ijo  is  large  compared  with  2/t, 


Vo+2/t  _ 

IC2-C-0I  ^  2  ' 


and  again  we  have  the  simple  bound 

l«-ol  s  l/HAioexp,)  -  Aigexp,!  <  «— Aiexp,. 

Vo 


By  continuing  this  process,  (3.9.5)  can  be  checked. 

When  the  q,  are  large  compared  with  n/r,  the  coefficient  of  1/2  that  appears  above  sug¬ 
gests  the  assymptotic  estimates  for  the  divided  differences,  as  in  Fig.  3.9.4. 

Estimated  error  bounds:  When  each  tj,,  7—0, 1 . //.is  large  compared  with  ii/t , 

|./f(Aii*'expr)  -  lij+'exprl  —  «K2~l(H-»i()_'Af0exPr  (3.9.6a) 

/-  o 

and 

l7/Uil+,expr)  -  Aii+/expr|  -  «r2-u+"(IIif<)-,-l/11+'-,eJtpr  (3.9.6b) 

,-n 

for  each  k  —0, 1 . //  and  /-2 . n—k  r\. 


(3.9.5)  shows  that  e«  times  a  matrix  like  that  in  Fig.  3.9.3  bounds  the  error.  (3.9.6a-b) 
indicates  that  ck  times  a  matrix  as  in  Fig.  3.9.4  is  a  good  approximate  bound.  The  elements  in 
Fig.  3.9.4  are  the  estimated  values  for  the  conjugate  differences  (3.7.14).  These  are  good  esti¬ 
mates  when  the  i\,  are  large;  hence,  a  bound  using  them  is  nearly  a  relative  error  bound. 
Because  k  depends  only  on  the  scaling  and  squaring,  the  standard  formula  portion  of  the  hybrid 
algorithm  does  not  lead  to  error  growth,  which  is  the  purpose  of  reordering  the  data  points. 

example.  The  data  from  the  previous  example,  namely 

i±0  -  -1.414214  ±  78.585786 

-  1.412799  ± /l  1.41563 
C±2  “  1-414214  ±  71 1.41421 

-  1.417039  ±  711.41 1 38  . 

generate  the  divided  differences  shown  in  Fig.  3.9.6  (only  half  the  table  is  exhi¬ 
bited).  The  data  were  generated  by  assigning  each  (,*°a  +  pe  /- 0, 1.2.3.  and 
rounding  to  seven  digits,  a  —  107,  p  — 2,  and  <6o“-3rr/4,  A|  —  tt/4  +  0.001,  <6i  —  rr/4 
and  <63 -ir/4- 0.002.  This  yields  both  closely  clustered  and  moderately  separated 
data  points.  The  arithmetic  is  seven  digit  single  precision,  so  condition  (3.4.14)  with 
spread  9  —  4  gives  7—2  squarings.  From  (3.9.2)  the  error  coefficient  is 

k  =  (3+1)(k,-4—  11  =  337. 


In  Fig.  3.9.5  the  estimated  bounds  are  very  close  to  the  true  absolute  values  of  the 


divided  differences,  so  k  indicates  a  loss  of,  at  most,  2.5  decimal  digits.  This  k 
clearly  excessive.  Indeed,  the  double  precision  coefficient 


k  =  *2-4-l  =  32.3. 


giving  a  loss  of  about  1.5  digits,  is  also  larger  than  the  results  in  Fig.  3.9.6  warrant. 


Differences  correct  to  seven  decimal  digits 

Row  Index 

l 

Divided  difference  table 

-1 

(  1.063 108E-3,  0.0  ) 

-0 

(  2.106638E-2,  0.0  )  (-5.2691 39E-2,  1.21 3605E-2) 

0 

(-1.624537E-1,  1.808715E-1)  (-3.7063 10E- 1,-1. 01 9594  ) 

1 

(  1.675059  ,-3.750361  ) 

-3 

(  1.61 1337E-6,  0.0  ) 

-2 

(  1.671 230E-4,  0.0  )  (  3.230809E-5,  1.838758E-5) 

-1 

(-1.1 14105E-4.  1.907577E-3)  (-2.524930E-4,  7.375033E-4)  j 

-0 

(-4.248672E-2.-2.540785E-3)  (-1.694748E-2.-3.852958E-3) 

0 

(-1.220463E-1.-8.447846E-1)  (-1  342108E-2.-3.523510E-1  >  1 

1 

(  1.673579  ,-3.754205  )  (  8.355560E- 1,-1. 880302  )  j 

2 

(  1.672096  ,-3.758050  )  (  1.669130  .-3.765728  ) 

3 

(  1.666154  ,-3.773412  ) 

-  ,  _  —  —  ....  J 

Differences  computed  by  hybrid  algorithm  j 

Row  Index 

Divided  difference  table  J 

-1 

(  1.063 106E-3,  0.0  ) 

-0 

(  2.106639E-2,  0.0  )  (-5.269139E-2.  1.213603E-2) 

0 

(-1.624537E-1.  1.808715E-1)  (-3.706307E-1.-1.019594  )  | 

l 

(  1.675059  ,-3.750361  1  1 

| 

-3 

(  1.61 1334E-6,  0.0  )  ! 

-2 

(  1.671228E-4,  0.0  )  (  3.230808E-5,  I.838754E-5) 

-1 

M.U4093E-4,  1.907575E-3)  (-2.524923E-4.  7.375030E-4)  j 

-0 

(-4.248669E-2.-2.540757E-3)  (-1.694747E-2.-3.852943E-3)  j 

0 

(-1.220468E-1.-8.447842E-1)  I-1.342132E-2.-3.523509E-1 ) 

1 

(  1.673576  ,-3.754205  )  <  8.355546E- 1 .- 1 .880302  ) 

2 

(  1.672096  ,-3.758050  )  (  1.669127  ,-3.765729  ) 

3 

(  1.666154  ,-3.773412  ) 

Fig.  3.9.6:  Conjugate  pair  exponential  divided  differences. 
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Exponential  divided  differences  with  real  abscissae  are  accurately  computed  by  a  hybrid 
type  algorithm.  This  idea  of  decomposing  the  divided  difference  table  into  blocks,  each  best 
computed  by  a  particular  method,  may  be  extended  to  additional  cases,  such  as  conjugate  pair 
exponential  divided  differences.  Indeed  any  sequence  of  abscissae  which  readily  decomposes 
into  well-separated  clusters  is  well  suited  to  the  hybrid  approach;  and  the  idea  need  not  be  res¬ 
tricted  to  exponential  divided  differences.  Though  scaling  and  squaring  does  not  work  in  gen¬ 
eral,  the  function  may  possess  special  properties  which  are  exploitable  through  representing  its 
divided  difference  table  as  a  matrix  function.  The  series  algorithms  are  still  applicable  for 
clustered  abscissae.  Certainly  many  extensions  are  possible,  only  the  simplest  and  most  basic 
have  been  dealt  with  here. 

Our  original  intention  in  studying  divided  differences  was  to  find  a  quick  and  accurate  way 
to  compute  the  matrix  exponential.  We  have  always  kept  in  mind  the  Newton  polynomial 
representation  and  techniques  appropriate  for  computing  matrix  functions.  The  techniques  we 
have  employed,  scaling  and  squaring,  the  standard  divided  difference  recurrence,  Taylor  series, 
and  decomposing  the  table  to  apply  a  hybrid  algorithm,  all  have  analogues  appropriate  for  com¬ 
puting  the  exponential  of  a  matrix  [Moler  and  Van  Loan,  I978I.  Indeed,  it  was  these  analogues 
that  suggested  many  of  the  approaches  pursued  here.  Thus  our  study  of  divided  differences  not 
only  aids  in  computing  more  general  matrix  functions  (the  Newton  polynomial),  but  it  also  pro¬ 
vides  an  indication  of  difficulties  that  lie  in  wait  in  matrix  function  evaluations.  Precisely 
because  divided  difference  tables  are  matrix  functions,  a  full  understanding  of  methods  for 
computing  such  tables  is  essential  to  an  understanding  of  functions  of  a  matrix. 
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A.l  The  Newton  divided  difference  series. 

For  those  readers  who  may  be  unfamiliar  with  divided  difference  expansions  such  as 

/<{>-  £  iovn  <{-<*,). 

k  — 0  /HJ 


we  present  here  a  convergence  proof  sufficient  for  our  purposes.  Similar  expansions  are  stu¬ 
died.  for  example,  by  Gel'fond  [1971],  but  they  are  not  quite  what  we  require. 

A  simple  derivation  of  the  Newton  divided  difference  series  is  obtained  from  the  contour 
integral  formula  (2.1.13) 


A0Y 


1  C _ f(to)  dot _ 

2iri'[.  (u  —  ao) (to  —  a i)  •  •  ■  (<o  —  a„) 


Our  proof  follows  a  method  commonly  employed  to  establish  the  convergence  of  complex  Tay¬ 
lor  series.  The  Taylor  expansion  of  /,  of  course,  is  a  special  case  of  the  more  general  Newton 
expansion. 

We  begin  by  deriving  a  Newton  formula  with  remainder.  The  expansion  points  are  the 
abscissae  of  the  divided  differences  which  are  coefficients  in  the  expansion. 

Newton  divided  difference  expansion  with  remainder.  Let  A„  =  {a0.al,  ...  ,aM)  be  a  se¬ 
quence  of  expansion  points  and  let  /  be  holomorphic  on  a  simply  connected  region  D  con¬ 
taining  A„.  Then  for  any  simple  closed  contour  C  in  D  enclosing  A„  and  a  point  £, 

/«)-  iAo‘/-ff  ({-«>>  +  «*({> 

*-0  ./- 0 

where  the  remainder 

*„({)  -  dto. 

ilrl  c  j-0 


«-t 

proof:  From  (2.7.5)  where  p„(0 

Jm  o 


p„+i((o)  - P„+M) 


1 


inu-op-nu-a/)) 

/-0  /- 0 


“  jinu-a,)-  n  (w— «()i , 

7-0  /-II  /-*+i 


Dividing  by  /Wi<m)  and  rearranging  yields 

“7  -  iin({-a,)-n(w_o',"'5 +  -«,><«-«, >- 

w  4  7-0  /-o  /-o  ®  fc  /-o 

By  Cauchy's  integral  formula, 

no  -  2^/(w  -  {)”'./ («)  du 

7-0  c  /-o  /-0 

-  i  A0‘/n  ({-«>)  +  *„<{).  □ 

7-0  7-0 

When  A„  consists  of  the  eigenvalues  of  a  matrix  -4, 


2flr/*£  (<u  —  4>(<u~ao)  •  •  •  («—  a„) 

is  the  characteristic  polynomial  of  .4.  When  /  is  holomorphic  inside  and  on 

/—a 

C,  the  integral  is  bounded  in  4.  R„(A)-0  by  the  Cayley-Hamilton  theorem,  thus  establishing 
the  Newton  polynomial  representation  of  f(A )  for  holomorphic  / 

We  need  only  show  that  /?„({)  —  0  as  n  —  oo  to  establish  the  Newton  series  formula. 


Newton  divided  difference  expansion.  Let  A  s  {ao,at|,a2,...}  be  a  sequence  of  expansion 
points  and  suppose  only  finitely  many  points  of  A  lie  outside  a  circle  of  radius  «  about  a 
point  a.  Suppose  further  that  /  is  holomorphic  on  a  simply  connected  region  D  containing 
A  and  a  disk  about  a  of  radius  p  >  2«.  Then  for  all  {  such  that  |{  -  a  |  <  p  -  2e, 

/(o  - 

k  — 0  /«0 


proof:  Select  a  simple  closed  contour  C  in  D  enclosing  A  and  such  that 
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Because  /  is  holomorphic  on  C.  there  exists  a  constant  K  such  that  |/(o»)|  <  K  for  ail  w  €  C. 

Let  M  s  \j  |  c  <  |a,  -  a |,  7-0. 1.2.... )  and  let  m  be  the  number  of  elements  in  M.  Since  m  is 
finite. 


0(f)  s  max  {1,  max  max  -p — ^4-| 

/€M  mil'  I Ctf'™  a  1 1 

exists.  C  was  selected  so  that  a*  ^  a,  for  any  j.  C  was  also  selected  such  that  for  all  m  €  C 

|«~{|  >  |w-a|  -  |f  —a)  >  2«, 
and 


I 


|a*-a,| 


p<— « 


=  y(f)  <  1 


for  all  j  €  M*.  the  complement  of  M.  Thus 


|/?„(f)|  < 


■\du>\ 


where  L  is  the  length  of  C.  Then  a s  /»  —  <»,  |/?„(C)|  —  0.  □ 

On  every  closed  disk  If  —  a|  ^p  where  p  <p  — 2e,  the  series  converges  uniformly  to  / 
p  may  be  chosen  as  the  radius  of  the  largest  open  disk  about  a  in  D.  When  the  sequence  of 
expansion  points  !«o.  «i.  «2.»1  converges  to  a,  the  «  of  the  theorem  may  be  chosen  arbitrarily 
small.  Convergence  of  the  Newton  expansion  may  then  be  claimed  for  all  f  such  that 
If  **  I  <P-  In  particular  when  all  the  expansion  points  are  confluent  at  a,  the  Taylor  expan¬ 
sion  of  /appears  as  a  corollary. 


) 

i 


Taylor  expansion.  Suppose  ./'is  holomorphic  on  a  simply  connected  region  D  that  contains 
a  disk  of  radius  p  about  a.  Then  for  all  (  such  that  |{  — a|  <  p. 


/<{>- 

*•0  *  • 


proof:  Recall  that  Ao/“. /■“’(<*)/*!  for  confluent  abscissae.  □ 

It  should  also  be  noted  that  because  /  is  holomorphic  on  D,  the  theorem  applies  equally 
well  to  any  derivative  of  /. 


,.2 


S 


4.2  Divided  difference  expansion  of  matrix  functions. 

The  results  of  the  previous  section  will  now  be  extended  to  functions  of  a  matrix.  For  a 
(w+1)  x  (//+!)  matrix  A ,  the  matrix  function  /(A)  has  a  Newton  series  representation  when 
,4’s  eigenvalues  all  lie  inside  the  series'  circle  of  convergence. 


Newton  divided  difference  series  for  a  matrix.  Suppose  ./'  has  a  Newton  series  expansion 
(as  in  §A.I )  on  the  disk  D,,  =  (£  |p-2«  >  |{  -a|).  Then  if  every  eigenvalue  4„  0  <  /  <  n, 
of  A  lies  in  D„. 


f(A)  -  £ AoV-II^ -«//>• 

A— 0  7*0 


proof:  For  any  at  ^  A(,  0  <  /  ^  /?,  the  matrix  («/  -  A )  is  non-singular  and 

(o»/-/4)_l  —  — a>/)-n(w  — a/)-1)  -I-  fj - —  (<o / -A)~' . 

*-0  ,-0  /- 0  0 


By  the  Cartan  definition  (1.1.7) 


f(A)  «  -y—f  -  A)~[d(o . 


The  simple  closed  contour  C  is  selected  such  that  it  encloses  ail  the  expansion  points  and 

pc  =  min  |w-o|  >  max  | A ,  —  at |  +  2e . 

w€C 


/(a )  -  ~aJn  +  r"(a) 

k  “0  7  *0 


where  the  remainder 


R„(A )  -  -r-ff(o>H<ul-A)-'n- 
2iriJc  fJo 


To  complete  the  proof,  we  need  only  show  that  in  some  norm  ||  R„(A )  |j  —»0  as  n  —  oo. 
Define  the  set  M  as  in  the  proof  in  §A.l.  Then1 


t  ||ff|L-  max 
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iirv  (u  to  — a.  ....  to  — a. 


/EM 

/<» 


/EM* 

/<« 


The  curve  C  has  finite  length  L  For  alt  to  €  C, 

I/M  <  K 

for  some  constant  K  because  /is  holomorphic  on  C,  and 

||(w/-/0-'IL  <  K' 


for  some  other  constant  K‘  because  C  is  bounded  away  from  A's  eigenvalues.  The  constant 

0  s  max  |1,  max  max  ||  — — —  IU} 

./EM  wEC  (O  — a, 

exists  because  M  is  a  finite  set  and  C  does  not  contain  any  a,.  For  all  j  €  Mc  each  eigenvalue 
X,  of  A,  0  <  /  <  «,  satisfies  the  inequalities 


lx,-«/l  ^  |X,-a|+<  „  _  |X,  —  a 

max -  <  -  ^  max  - 

«ec  to  — a/ 


+  e  , 

max  — - - - =  y  <  1 . 

Pc“*  0<i*»  pc- 4 


Let  A  -  P~'JP  where  the  upper  triangular  matrix  J  is  /4’s  Jordan  canonical  form.  Then 

n  -  p-'-  n 

*  £U  —  Of  .  **  O)  —  Ct , 

/EM*  >  /EM*  1 

I  <  »  /  <  » 


-  p-'d •  n 


D~'JD-«,l 


j  EM* 

/<» 


tO  —  Ctj 


D~'P 


where  £>  =  diag(l.T},T>2. . .  .  ij  >  0.  Taking  norms. 

ii  n  ^  ii/’-,z>iu-i£>-,/,iL-  n  ii 

(tl  —  Q  ■  6)  *  Q  , 

JEM*  '  /EM*  J 

y<» 


EUlZll ||.  <  m„  +  -JL 

to  — at)  o </<«  |»— a;  way 


<  y  +  — *  y  <  1 

Pc“« 


For  any  j  €  Me 


for  all  -»i  <  (pt — «)<1  — y')-  Thus 


ii  n  ^^ii-  <  *v- . 


where  the  constant  K"  =  ||/>-lD||«-||Z>-,/,||»  for  some  fixed  tj  <  (pc  —  « ) ( 1  - 
these  bounds  yields 


and  || /?„(/<) I!,.— 0  as  n  —  □ 


’).  Combining 


1 
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B.  Error  bounds  for  a  Taylor  series  computation  of  4<?expr. 

The  Taylor  series  is  one  of  a  number  of  proposed  methods  for  computing  4n'expr.  It 
involves  no  previously  computed  divided  differences;  so  its  study  does  not  concern  propagation 
of  errors,  but  just  computational  errors.  Here  we  develop  error  bounds  on  the  computation  of 
Ao'expr  by  a  Taylor  series  and  demonstrate  that  the  method  is  best  applied  when  the  abscissae 
are  closely  clustered. 

In  §2.8  A<j'expr,  with  real  or  complex  abscissae  {{o.£i.  ■  •  •  .{„)  and  r  >0,  is  shown  to 
have  a  Taylor  expansion  about  a 

StSSt*""'-  (BU 

where  the  power  function  f "+/  is  t„+/({)  -  (£  -  a)"+',  j  -*0, 1 . 2 .  It  is  convenient  to  consider 

the  shifted  abscissae  ({0-a,{|-a . £„-«)  exact;  the  numerical  effects  of  shifting  abscissae 

are  discussed  in  §3.6.  With 

S  =  max  |£,-a| .  (B.2) 

0<i<» 

the  bounds  we  obtain  resemble 


|/f(i<?expT)  “  4o'expr|  <  tie 


where  At  represents  a  coefficient  dependent  on  the  arithmetic  details  to  be  introduced  shortly. 
V?(4d'expr)  represents  the  computed  floating-point  value  of  Ao'exPr- 

The  Taylor  series  algorithm  outlined  in  §3.3  requires  many  inner  products.  We  consider 
two  separate  conditions  for  bounding  round-off  error  accumulation  in  inner  product  computa- 
tions.f 

M 

I.  Double  precision  accumulation.  The  error  in  the  computed  inner  product 

!m 0 

satisfies 


<  cfk/3,1 .  (B.3) 

f—0  t— 0  0 

tSee  Wilkinson  (1963)  for  a  general  treatment  of  rounding  error  analysis. 
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where  «  =  1.06  x  (machine  precision),  c  is  assumed  so  small  that  any  0(<2)  expressions  are 
negligible  when  compared  with  expressions  linear  in  <  and  may,  because  or  the  arbitrary  1.06, 
be  absorbed  into  such  linear  expressions.  Error  condition  (B.3)  holds  for  double  precision 
accumulation  of  sums  and  inner  products.  It  does  not  depend  on  the  number  of  terms  summed 
and  leads  to  simple  illustrative  error  bounds.  Additionally,  we  assume  the  series  coefficients 
r"*'eTa/(n+j)[  are  all  calculable  to  machine  precision,  namely 


[fill 


H+J.TU 


TM+ieTa  J  _  T'^'e 


<  c 


r"T'  \e' 


( n+j)\  ( «+/)!  («+/)! 


(B.4) 


2.  Single  precision  accumulation.  The  second  condition  applies  to  single  precision  computation 
of  all  quantities.  Wilkinson  [1963]  shows  that 


L ~  <  «[U+l)|ao0ol  +  £(«+2-/)|a,/3,|] . 

0  /— 0  /—I 


We  simplify  this  to  the  more  convenient 


“  ±afi,\  <  «£(»+2-/)k/8,|  •  (B.5) 

f“fl  i 


In  addition,  we  assume  the  series  coefficients  are  evaluable  with  no  more  than  five  rounding 
errors  (say,  errors  in  the  evaluation  of  r"+/,  eTa,  and  (n+j)\,  plus  a  multiplication  and  a  divi¬ 
sion);  so 


I 

•f 


i 


u 


(B.6) 

Bounds  derived  from  (B.5),  though  more  complex,  are  more  generally  applicable  than  those 
from  the  first  condition. 


We  start  by  deriving  bounds  on  divided  differences  of  power  functions. 


Lemma  1 :  For  y  — 1,2,... 

.  ,  n . 

(B.7) 

proof:  From  the  recurrence  (2.7.8), 


-  I({,-a)-16tr'''- 

i-0 


For 7-1.  <fcol«+l"S(t'-Q[*  and  s° 

i—O 


Uoti+li  <  (k+\)&.  k- 0,1 . /;. 


If  for  some  j  >  1 


for  each  k  -0. 1 . n,  then  since  Aoti+t-£(£,-a)-A()t<'+'~1, 

,-o 


IAa‘T,+*1  <  fi7  *r  (/+/-1)’  (J+k)  \ 

!A°L  1  *  (7-1)!  J  /!  kTJl  6 


for  *-0,1 . n.  □ 


We  now  give  bounds  on  the  error  in  jKlfta**)  f°r  each  k  and  7.  When  the  error  is  not 
too  large,  7?(Aoi„+'‘)  may  also  be  bounded  as  in  Lemma  1. 


Lemma  2:  Let  AoU+*  be  computed  according  to  Algorithm  1  of  §2.7.  Then  for  each 


7  — 1.2,...  and  *—  0, 1,  .  .  .  ,n. 


W*4\i+k)  -  Ao‘T«i,+‘l  <  -  j**J  J\k 


w-C'-uio  '  -uia  1  ^  k\(j-\)  I  j  «  ) 

for  double  precision  accumulation  (B.3),  while  for  single  precision  accumulation  (B.5), 


(B.8a) 


l/WAo*ir*>  -  A0‘irAl  <  +  o-n-^7 ~-i  • 


(B.8b) 


proof:  For  (B. 8a), 


lyWta‘+'>-A<frt«‘+,l  <  «£|{,-«l  <  <*+!>««.  *-0.1. 

1-0 


by  (B.3)  when  7  —  1.  If  for  some  j  ^  1 


Theorem:  In  computing  Ao'expr  by  a  Taylor  expansion  about  a,  the  error  is  bounded  by 


L/lj(Ao'expr)  -  A0"expT|  <  cer*(2-HrS)  T  I1 ■ 

n ! 

for  double  precision  accumulation  (B.3),  while  for  single  precision  accumulation  (B.S) 


l//(Ao'expr)  -  Ao'expT|  <  «(m  +  n  +  7  +  rS)er*-^l^- . 

n\ 

where  w+ 1  is  the  number  of  terms  actually  summed  in  the  expansion. 


(B.9b) 


proof:  Let  A(j' s,l+nl  =  ]jT/3,l+>-Aot «+'  be  partial  sums  of  the  Taylor  expansion  (B.l),  where  each 

/- o 

coefficient  /3„+(  =  T"*'era/(n+j)\.  The  error  is  bounded  by  four  terms. 


L//(Ao'expr)  -  Ao'expr|  <  lftUfexpT)  -  £/W0,I+,)./KA(j' |"+')( 


+  IL/K/W^Aoir'>  -  I/3„+^(Ao'l''+')| 
/— 0  /“0 


+  H/3„+,-yKA<nr')  - 

Jm 0  ;—0 


+  |Ao's„+m  -  Ao'expr| 


=  I  +  II  +  III  +  IV. 


We  bound  each  of  the  four  terms  separately.  In  addition  we  note  that 


/f(Ao'expr) 

/•o 


Double  precision  accumulation  (B.3):  By  (B.3)  and  (B.7), 


i  «£i/2</wiL/MA<Ftry>i  < 

jm o  ;-0  n  J 


n\  j\  n\ 


By  (B.4)  and  (B.7), 
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1  \ 


/-o 


1-0 

And  using  (B.8a) , 


»  .w|-,ra| 


n\ 


in «  - «tn  < «ir(* I')’1 


jUL.  s 


We  may  ignore  the  truncation  error  IV  because 

IH  <» 

iv  <  n/w^tr'-  i/wa*?:+/i 


/-o 


y-0 


<  £  i/s^jiAd'Tri  <  1,7 


/•w+l 


/  —  Wf  +  I 


is  negligible  for  m  large  enough.  Summing  the  bounds,  then,  yields  (B.9a). 

Single  precision  accumulation  (B.5):  The  same  steps  are  repeated  for  (B.9b).  By  (B.5)  and 
(B.7), 


1  <  t Z (m+2-y)  l/7(j8w+/)  |  l/WAd’T 2+')  |  <  «f ,  (w+2-y) 


/-0 


J-  0 


rW 


In  the  same  manner  as  before,  by  (B.6) 


II  *5  5«<?T*-^ 


"l-j’Oj 
I 


n\ 


By  (B.8b) 


III  s$ 


-  utr^i 


i-0 


+  Q-D  frtS)L}y 

*1  <«+»!  y»!0+D!  V  n!y!  ' 


U 
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■  t- 


and  so 


,  ,  ^  T"\era\  ,  ,  ,  jn  ,  r'V 


<  eitn  +  n  +  1  +  r8)er*- 


Finally,  we  choose  in  so  large  that,  say, 


£ 

y-«+ 1  J  • 


(B.10) 


From  the  discussion  earlier  Tor  IV, 


n\ 


Summing  our  bounds  on  I  +  lII,  II  and  IV  yields  (B.9b).  □ 


r"eTa/nl  is  Ao'expr  for  abscissae  confluent  at  a,  the  bounds  (B.9a-b)  make  clear  why  the 
Taylor  series  method  is  best  applied  to  closely  clustered  abscissae. 

Relation  (B.10)  permits  determination  of  m  when  a  particular  c  and  t8  are  given.  For 
example  when  «- 10-7  and  r8  <  l,  (B.10)  yields 

£  4r  ee1  =  2.72  xlO~7. 


The  smallest  value  of  m  for  which  this  inequality  holds  is  m  - 10.  And  when  «  —  10~IJ.  the 
smallest  m  is  m  —  1 6. 


C.  Decision  criteria  for  the  hybrid  algorithm. 

Double  precision  accumulation.  In  §3.5  we  found  that  the  decision  criterion  rt»„  for  the  hybrid 
divided  difference  algorithm,  with  double  precision  accumulation  (3.3.3),  satisfies  the 
recurrence  (3.5.2) 


(k2t9„-  I)  -  (kjt#„_|-  1 ) ( 1  +  2«/r0„)  -  0  (C.l) 

where  r«„  =  2/k2.  This  recurrence  has  no  simple  closed  form  solution  for  r«„.  However  it  is 
possible  to  give  a  simple  bound  on  t9„. 

The  recurrence  (C.  1 )  is  quadratic  in  r#„;  thus 

r«„  -  y{T0„_,  +  >/(t9h-,)2  +  8n(r«„.l-l/»t2))  (C.2) 

is  a  rearrangement  of  (C.l)  where  t9„  appears  explicitly.  We  attempt  to  bound  t9„  for  every  n 
by  linding  a  function  in  n  which  satisfies  a  majorizing  recurrence.  A  little  exercise  in  complet¬ 
ing  the  square  gives 

/i(/f+1)+2/kj«‘‘  y{»(»—l  )+2/k2  +  ^[rt  («~l)+2/*c2l2  +  8/i(/r(n— 1)+2/k2— l/itj]  +  16nJ  +  8«/»c2j  . 

which  is  nearly  the  same  as  (C.2).  Since  t0o-2/k2,  it  is  dear  that 

t9„  <  n(n+l)  +  2/K2  (C.3) 

for  all  n  >  0  and  any  k2  >  0.  Also  in  a  similar  way, 

it (n-. 3)  -  yf(*-l)(*-4)  +  y/l(n-\)(n-4)}2  +  8nUn-\)(n-4)  -  \/k2]  +  16/7  +  8«/*2}. 

We  compare  this  with  recurrence  (C.2).  For  the  value  of  k2- 8.3259  derived  in  §3.4.  we  find 
that  r0|7- 237.85  <  17  (17 -3) » 238  from  the  table  in  Fig.  3.5.3.  Thus  t0„  <  n(n- 3)  for  all 
n  >  17.  However, 

«(fl— 4)  *  y  {(«— 1 )(«— 5)  +  l)(n-5)|J  +  8/t  l(rr— 1)(«— 5)  —  1/k2]  —  4«2  +  40n  +  8n/*c2) . 

For  «2-8.3259,  -4ff2  +  40ff +  8«/k2  <  0  when  n  >11.  Since  r0io- 72.02  >  KM10-4) -60. 
t9„  >  n(n- 4)  for  all  n  >  10.  Combining  these  two  results  yields,  for  n  >  17  and  k2-8.3259, 
that  t9„  is  bracketed  by 


— --  -  a m 


n(n— 4)  <  t9„  <  n(n—  3). 


(C.4) 


Single  precision  accumulation.  The  decision  criterion  rtf*  for  single  precision  accumulation 
(3.3.4)  satisfies  the  recurrence  (3.S.6),  namely 

(/t+!)(K,r0H-l)  -  —  1)(1  +2n/r9„)  —  0.  (C.5) 

This  is  also  quadratic  in  r«„;  so  we  have  the  equivalent  recurrence 

r»„  ■»  2k  (n+  jj U* i » Tg„-i  +  1  >  +  >/ (K|/rT©„_|  + 1)2  +  8K|rtJ(w+l)(K|T0M_|  -  l)|  (C.6) 

in  which  r0„  appears  explicitly.  Initially  t#i -3/2*!. 

For  cr„  =  2/»J/3  +  (3/2k|  -2/3)»,  we  find  by  completing  the  square  that 

<T"  “  2k  (1+1)  +  +  8k,/i2(/H-1)(k|0-„.|-1)  +  m,,)  (C.7) 

where  ff„_|-2(/i-l)J/3  +  (3/2K|-2/3)(n-l)  and 

v„  —  4k|(/t+1){(2  +  +  (-r- - r—  +  -~)/j2  +  (-i-  —  -r — )»} . 

3  4k  §  9  3  3  2k  j 

o-„  was  chosen  so  that  <n  —  3/2ki  —  t0|  and  v„  >  0  for  any  K|  >  0  when  n  >  2.  Comparing  the 
recurrences  for  <r„  and  r0„  shows  that 

r0„  <  ~n2+  i^-hn  (C.8) 

3  2k|  3 

for  n  >  1. 

To  bracket  t9,„  for  large  n,  <r„  =  n(2n-S)/3  satisfies  the  recurrence  (C.7)  with 
V„  —  4K|(ff  +  l){(2Kt  +  —  )/t2  +  yfl) . 

A  check  of  the  first  and  fourth  columns  of  Fig.  3.5.4,  which  has  k,- 21.2950,  reveals  that 
3.16- t04  <  (r4-4.  Thus  t9„  <  <r„  SB  a(2/»-5)/3  for  all  n>4.  Similarly,  <r„  3  2n(n-3)/3 
satisfies  the  recurrence  (C.7)  with 

v„  <■  4k)(«+1)|-- ~ -n3  +  (4k|  +  y )n2  +  2n | . 
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K 

When  n  >  10,  v„  <  0.  From  Fig.  3.5.4,  120. 12  — r9,s  >  o-,,- 120.  Thus  >  2«(«-3>/3  for 
n  >  15.  Combining  the  two  bounds  shows  that  rtf*  is  bracketed  by 

jn2-2n  <  t«h  <  (C.9) 

for  n  >  15  and  k,  -21.2950. 
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D.  Numerical  examples. 

The  tables  on  the  following  pages  illustrate  the  example  in  §3.6  of  the  hybrid  algorithm 
with  clustering.  The  first  table  (3  pages)  is  the  hybrid  compulation  in  single  precision  for  r  - 1 . 
The  correct  seven  digit  divided  differences  are  presented  in  the  following  table  for  comparison. 
The  two  following  tables  exhibit  in  a  digits  lost  (logi0>  form  the  actual  relative  error  and  the 
results  of  an  a  priori  error  bound  computation.  The  data  in  these  tables  are  summarized  by  Fig. 
3.6.3.  A  second  set  of  tables  for  r  -  2  then  follows  (see  Fig.  3.6.4).  Finally  for  comparison, 
the  table  for  r-2  is  recomputed  by  scaling  and  squaring  only  (Fig.  3.6.5).  The  abscissae  are 
listed  to  the  left  of  each  table.  The  computations  were  performed  on  a  PDP-11  computer, 
which  has  a  precision  slightly  greater  than  seven  decimal  digits. 
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